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mathematical concepts and historical development of the 

MASCOT GUIDANCE TECHNIQUE FOR SPACE VEHICLES 

SUMMARY 


The recent development of improved integration routines, improved 
trajectory computation algorithms, and faster, lighter, more flexible 
flight -worthy digital computers has opened up new possibilities for 
improved guidance concepts and effective application of optimal guidance 
schemes to meet the challenges of the post-Apollo space vehicles . 

This report documents one such guidance concept which permits 
practical realization of many heretofore unrealizable optimal guidance 
tasks, namely, on-board optimization of ascent trajectories through the 
atmosphere, self-targeting, optimal rendezvous, optimal coast-burn- 
coast-burn orbital transfer, and reentry for high lif t-to-drag-ratio 
type of vehicles. 

This collection of techniques combined into one guidance program 
has been called MASCOT (Manned Shuttle Comprehensive Optimization and 
Targeting). 


I. INTRODUCTION 


The advent of a new generation of space vehicles (specifically 
NASA’s post-Apollo space shuttle) with advanced features such as com- 
pletely reusable stages, sophisticated on-board digital computers, and 
airplane- type landing characteristics for both booster and orbiter has 
required a re-evaluation of the space vehicle guidance concepts used 
in the past. When re-examined, these traditional guidance concepts 
have been found to be inadequate for the new generation of space 
vehicles. In response to this need, new concepts, procedures, and 
techniques have recently been developed which provide a significant 
improvement in optimal guidance applications. In particular, this new 
approach yields a physically realizable optimal guidance scheme which 
results in nearly optimal performance for the boost phase, the rendezvous 
phase, and the reentry phase of the post-Apollo shuttle vehicle flight 
path. This new system and concept has been given the name MASCOT. 



The MASCOT concept is made possible by the development of new 
mathematical techniques and the significant improvement of others. It 
also takes advantage of the improved speed, memory capacity, and relia- 
bility of the new flight-worthy digital computers * 

In this report, the evolution and mathematical theory of the MASCOT 
guidance technique are explained. In addition, some computer-s imulated 
trajectories obtained by means of the MASCOT program are presented to 
demonstrate the accuracy, flexibility, and optimality of the MASCOT 
system. Some comparisons will be made, whenever practical, with 
results obtained from the present Saturn Iterative Guidance Mode (IGM). 


II. SOME BACKGROUND ON GUIDANCE, NAVIGATION, AND CONTROL 


The operation of a guidance system involves three interrelated 
divisions of effort, navigation, guidance signal generation, and con- 
trol, Navigation is concerned with the determination of the current 
position and velocity of the vehicle. Navigation is normally an open- 
loop process for short flights but must be made closed-loop by naviga- 
tion update in some instances, i.e., long flight times. 

With an accurate knowledge of all physical parameters of vehicle 
and environment, current position and velocity, and an accurate know- 
ledge of the desired terminal conditions, it is possible to determine 
the proper directions of the engines thrust to obtain the desired end 
conditions. This closed-loop process of thrust-direction determination 
is referred to as guidance signal generation. 

The control system implements the guidance signals in addition to 
employing such maneuvers that may be required to preserve the structural 
integrity of the vehicle. Thus, the control system may temporarily 
ignore or minimize the influence of the guidance signal in order to 
protect the vehicle from winds, wind gusts, wind shears, vibrational 
disturbances, or fuel sloshing problems. Control, too, is a closed- 
loop process. 


III. THE GENERAL CONCEPT OF GUIDANCE 


Implementation of the guidance signal will determine the vehicle 
flight path to reach a desired destination or set of terminal condi- 
tions. Perfect operation of the navigation system, exact generation of 
the guidance signal, and exact implementation of the control system, 
together with an exact knowledge of the vehicle physical parameters and 
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an exact knowledge of the physical environment daring the entire flight 
to the desired end conditions, would make it possible to solve for the 
guidance signal only once. This signal would, of course, be a function 
of time or some other convenient variable so that the guidance signal 
would change continuously during the flight. 

In practice, perfect operation of the three systems cannot be 
achieved. Thus, a practical guidance system may make the computations 
at the initiation of guidance, and predict the required guidance signal 
or function to guide to the required end conditions. Since the initial 
information will have some inaccuracies, the process must be repeated 
periodically during the flight, treating each instant at which the new 
computations occur as a new set of initial conditions. 


IV . OPTIMAL GUIDANCE 


Optimal guidance, a subclass of the broad class of guidance con- 
cepts, is characterized usually by the deliberate maximization or mini- 
misation of some aspect of vehicle performance, for example, maximum 
payload to a given orbit or equivalently minimum fuel required for a 
given payload. 

To accomplish optimal guidance of a vehicle, a mathematical model 
of the vehicle motions must first be chosen. Appearing in these equa- 
tions of motion are the so-called control variables. These variables 
are the equation parameters which can be given an arbitrary (within 
prescribed limits that depend on the problem) value at any instant along 
the vehicle* s flight path. This large amount of freedom of selection 
for the control variables allows a time history for the control variables 
to be chosen which achieves a desired destination and also optimizes 
some specified aspect of the vehicle performance. For a typical rocket- 
powered vehicle, the control variable is usually the thrust direction, 
and a time history of this thrust direction is determined (as described 
previously) to maximize the payload delivered to a particular destination. 


A. Mathematical Techniques of Optimal Guidance 

The different mathematical techniques used to select the optimum 
control variable are usually grouped into a body of knowledge known as 
optimization theory or calculus-of-variations theory. 

When some type of optimization theory is applied to the differen- 
tial equations simulating the motion of a rocket-powered vehicle, a two- 
point boundary value problem is produced. This means that not all of 
the initial conditions for the differential equations involved are known 
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at the initial time. The unknown initial conditions must be determined 
so that they cause the results of the integration of the differential 
equations to satisfy some pre-specif ied conditions at the final time 
or destination. To accomplish this, in practice, one usually attempts 
to discover, numerically, some mathematical relations which connect 
the initial conditions with the final conditions. Then, the desired 
conditions at the destination are placed in these relations, and the 
resulting equations are solved, hopefully, to yield values for the 
missing initial conditions. These initial conditions and the relations 
connecting the initial conditions with the final conditions yield the 
time history of the control variables that optimize the chosen criterion 
of the vehicle's performance. In practice, there may be additional 
relations or constraints to be satisfied by a particular vehicle's 
trajectory which produce a multi-point boundary value problem instead 
of the simple two-point boundary value problem described above. These 
aspects of optimization theory are examined in more detail in the sec- 
tion on mathematical development. 


B. The Philosophy of Optimal Guidance: Old and New 

In the past, computer technology and mathematical techniques were 
inadequate to allow the consideration of a realistic computational 
algorithm for trajectory optimization as a practical on-board guidance 
signal generator. Thus, older techniques for on-board guidance signal 
generation relaxed the reality of the mathematical simulation of the 
vehicle's motion or degraded the optimality of the solution path in 
order to solve the resultant boundary value problem rapidly enough to 
control a particular vehicle's flight. The new MASCOT concept (which 
is a realistic computational algorithm for trajectory optimization) can 
be considered as a guidance signal generator for all phases of a 
vehicle's flight (boost phase, rendezvous phase, and reentry phase). 
This advance is possible at the present time because of improvements 
in computer technology and mathematical techniques for the solution of 
optimal trajectory problems. 

Before going directly into the explanation and results associated 
with MASCOT, we briefly outline the historical developments which led 
to the present form of the MASCOT guidance technique. For readers not 
familiar with some of the ideas and terminology contained in this 
historical development, an extensive list of references andia biblio- 
graphy * are ihcludjedrat the end of the report. 
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V, HISTORICAL BACKGROUND 


A, Delta Minimum Guidance 

Beginning with the V-2 rockets in the mid-I940 f s and continuing 
through the Redstone, Jupiter, and Pershing missiles, the "delta- 
minimum" guidance concept was used successfully to guide rocket flights. 
The delta-minimum concept required that both nominal and perturbed tra- 
jectories follow essentially the same geometrical trajectory regardless 
of other considerations. On-board flight computations were simple, 
and analog computation was used for the execution of the delta-minimum 
concept. 


B. Iterative Guidance Mode 


In 1960, research work was begun to develop new guidance concepts 
for the Saturn space vehicles [1,2,3]. This work was motivated by the 
development of new mathematical techniques for maximization of payload 
through optimization methods and by the development of digital computers 
to replace analog computers as on-board hardware. It was also obvious 
that space trajectories would require greater flexibility to cope with 
sudden changes, such as engine out, and that more options must be per- 
mitted in the selection of flight profiles. The guidance law which 
resulted from the research work begun in 1960 and which has been flown 
on Saturn vehicles was given the name IGM, Iterative Guidance Mode. 

IGM is essentially an approximate formulation of the calculus-of- 
variations problem that allows analytic construction of major parts of 
the solution, so that only a simple iterative numerical process is 
required for solution. This approach avoids the time-consuming numeri- 
cal integration procedures that have been required to compute a general 
solution to the fundamental calculus -of -variations problems. The speed 
needed for real-time application has been the primary motivation for 
the derivation of semi-explicit methods of this type. 

However, as a result of the approximations, the optimality and flexi- 
bility of such flight schemes are limited, primarily in that they are 
nearly optimal only for short arcs of powered flights and for specialized 
mission (boundary value) conditions in a restricted coordinate system. 

This limitation can be relaxed somewhat in practice by use of special 
purpose adjustments, but only at the expense of additional preflight 
analysis and simulation. 

The basic simplifications made in the Iterative Guidance Mode to 
obtain analytic construction of major parts of the solution were (1) to 


5 



assume a uniform gravitational field rather than an inverse square law 
and (2) to apply the same transversality condition at engine cutoff 
regardless of the mission flown. There were, incidentally, some small- 
angle approximations and ingenious special purpose adjustments made to 
improve the optimality of the scheme and to reduce the computational 
effort required for on-board mechanization. However, the basic concept 
of the simplifications is best seen in view of the assumption concern- 
ing the gravitational field and the assumption of the adequacy of a 
single transversality condition. 


C. QUOTA 

Next came the QUas i Optimal Trajectory Analysis, QUOTA [4], which 
belongs in the same category with IGM since it too is an approximation 
to the calculus-of-variations (COV) solution. 

The fundamental approximation in QUOTA may be viewed in at least 
two different ways. The first way may be said to be an approximation 
to the COV in the sense that the Euler-Lagrange equations are replaced 
by a slightly different set of equations, the purpose of which is to 
allow analytic integration of both the "state" and "co-state" equations. 
Of course, the "co-state" variables must be interpreted in a generalized 
sense because they are not obtained from the Euler-Lagrange equations. 

A second manner of viewing the approximations to achieve the solu- 
tions will be mentioned. From many observations of the behavior of the 
Lagrange multipliers (Vs) in the calculus of variations, it has been 
noted that for rocket flights in a wide variety of missions, these A* s 
are very nearly linear functions of time. The QUOTA equations may be 
derived by assuming linear A f s and by expanding the gravity term in 
s eries . 

The advantage of this approach is a very rapid computational 
scheme which is much more flexible than IGM and also leads to a smaller 
loss of optimality. Even this small loss of optimality may be regained 
by making the assumption that the A f s are linear only for a portion of 
the flight; thus, the A versus time curve is represented by a series 
of connecting straight lines where the slopes of the linear A segments 
are determined from the Euler-Lagrange equations. For an Apollo- type 
mission, the A 1 s may be assumed to be linear for the entire flight with 
only a negligible loss in payload. 

With either of the viewpoints mentioned above, the true gravita- 
tional field may be accurately represented and transversality conditions 
may be fitted to the particular end conditions. Needless to say, the 
optimality of the solution must be checked because of the substitution 
of approximate equations for the Euler-Lagrange equations. 
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Computationally, QUOTA offered much more flexibility than IGM, 
required less preflight analysis and only slightly increased the 
amount of on-board computation required. Because IGM was entirely 
adequate for all Saturn flights, there has been no need to change 
from IGM to QUOTA and, therefore, QUOTA has never been tested on any 
actual flights. It is still available for any mission which requires 
flexibility beyond IGM capability and for which the same computational 
equipment must be used. 


D, OPGUID AND SWITCH 

Before 1965, general (flexible) numerical procedures for computing 
precise optimal trajectories were too unreliable in convergence and 
costly in computational requirements to be considered for real-time 
guidance. However, an indirect method for computing optimal trajec- 
tories [5,6] was developed in 1965 incorporating improved techniques 
to obtain a substantial gain in speed, convergence, and flexibility. 

A simple scaling rule for the amount of the Newton correction that was 
permitted per iteration resulted in an extremely large region of con- 
vergence that was surprisingly insensitive to accurate initialization 
of the boundary value problem. 

The indirect approach is particularly well suited for real-time 
use, where the continual adjustment of the guidance scheme to accommo- 
date perturbations in initial conditions is readily accomplished by a 
single Newton iteration on the boundary value problem. In 1966, the 
feasibility of this approach as a real-time guidance scheme for optimiz- 
ing single-burn-arc orbital injection missions was demonstrated [7] and 
named OPGUID (OPtimal GUID ance) . 

However, many orbit transfer problems require the use of several 
burn arcs separated by relatively long optimal coast arcs. A multi- 
burn-arc version of OPGUID, developed in 1967 [8], demonstrated that 
the attractive fundamental approach of OPGUID could successfully con- 
verge a general formulation of this problem, with variable boundary 
conditions. A sophisticated version of the multi-burn program (SWITCH) 
has been developed [9] that has successfully converged a variety of 
orbital transfer problems with an efficiency and reliability comparable 
to that of the original OPGUID. 

As a result, the indirect method of SWITCH is not only feasible but 
considerably superior to existing implementations of quasilinearization 
in convergence as well as efficiency. A principal feature of SWITCH is 
the use of classical two-body theory to render the computations for 
coast arcs explicit. Since high-thrust multi-burn orbit transfers'" 
usually involve coast arcs many times longer in duration than burn arcs. 
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the explicit method results in a substantial savings in computation 
per iteration. A universal variable formulation of the two-body 
problem with closed-form expressions for the state transition matrix 
is used. This formulation was adapted from the work in references 
8, 10, 11 and 12 in a novel way to avoid the cumbersome computation 
of the three-dimensional tensor of second partial derivatives of final 
state with respect to initial state that is required when computing 
the partial derivative of final co-state with respect to initial state. 

Since all the partial derivatives are available at each iteration, 
as was the case for the original OPGUID, the SWITCH algorithm is appro- 
priate for computing real-time corrections to in-flight perturbations . 

Unlike the OPGUID algorithm, the SWITCH algorithm does require 
reasonable initialization. That is, it is not possible with SWITCH 
as it was with OPGUID to misalign the thrust direction by 90 or 180 
degrees and retain convergence. However, rough estimates of impulsive 
solutions have proved more than adequate as initialization in every 
trial case. 


E. MASCOT 


With the advent of the space shuttle with its high aerodynamic 
lifting characteristics both on ascent and descent, it has become 
imperative that optimal guidance laws be devised to optimally steer 
the vehicle through the atmospheric portions of flight. 

This has been accomplished essentially by introducing into the 
SWITCH program the aerodynamic effects of both lift and drag [13], 
Additionally, the effect on thrust of pressure variation with altitude 
has been included in the propulsion computations. 

These changes have recently been carried out, and some of the 
observed results may be summarized as follows: 

1. Maximum payload and maximum sensitivity to convergence occur 
when atmospheric effects are introduced into both the state and co-state 
equations . 

2. Approximately two- thirds of the payload increases obtained 
under (1) may be kept by introducing the atmospheric effects into the 
state equations only. In this case, convergence is far easier. 

3. Extremely difficult convergence Cases may be approached 
gradually by multiplying the atmospheric model by 0 ^ R ^ 1 where K 
starts at zero and gradually increases to 1. For each value of K, the 
case is converged before increasing the value of K. 
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When the atmospheric effects, together with changes to be described 
later, were added to the SWITCH concept, the name MASCOT was chosen for 
the overall scheme to handle multiple burn-coast-burn optimal trajec- 
tories both in and out of the atmosphere. 


VI. MATHEMATICAL DEVELOPMENT FOR MASCOT 


A. The Choice of a Physical Model Common to All Flight Phases 

The first step in the development of an optimal guidance technique 
for a rocket-powered vehicle is the mathematical modeling of the vehicle 1 s 
trajectory. To do this, a second order three-dimensional vector of 
ordinary differential equations is usually used which consists of the 
sums of the vector accelerations produced by all the forces to be con- 
sidered acting on the vehicle. For the MASCOT scheme, this sum of vector 
accelerations must consist of terms produced by forces acting during the 
boost phase, the rendezvous phase, and the reentry phase. During the 
computation of a trajectory^ for a particular phase, all of the terms in 
the sum for all of the phases do not have to be considered and the 
unneeded terms can be easily ignored with the logic of the computer 
program. The next section develops in detail the general mathematical 
model, and then the terms to be ignored in particular phases are pointed 
out specifically. 


B. Mathematical Modeling of Vehicle 1 s Motion 
1. Comments on Coordinate Transformation 

A three-dimensional Cartesian coordinate system with an origin at 
the earth* s center will be used. The orientation of this coordinate 
system will be assumed to remain fixed with respect to time and the 
initial orientation is completely arbitrary. To compare trajectories 
in this arbitrarily oriented coordinate system with any other coordinate 
system, their relative orientations at comparison times must be known so 
that a coordinate transformation can be performed on one of the coordin- 
ate systems. Some coordinate systems of interest do not always have 
their origin at the earth* s center, but this is a simple translation if 
relative orientations are the same. For example, in booster work, an 
inertial coordinate system with origin at the earth* s center could be 
defined to have the z-axis perpendicular to the earth* s equatorial plane 
and pointing northward, the x-axis in the equatorial plane and pointing 
toward the Greenwich Meridian at some reference time (usually launch time), 
and the y-axis pjbaced in the equatorial plane to form a right-hand system. 
Then, any other coordinate system can be oriented with respect to this 
coordinate system with the aid of the standard Euler angle rotation 
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matrices. That is, for a given angle a , the Euler angle rotation matrices 
for rotations about the x, y, and z axis are written as follows: 



i 

0 


0 


II 

H 

Q 

0 

cos 

a 

sin 

a 


0 

-sin 

a 

cos 

a 


(a) 2 


cos a o 
o l 
sin a 0 


-sin Q£ 
0 

cos a 


(a) 3 


cos a sin a 0 

-sin a cos a 0 

0 0 1 


Then coordinates in the equatorial coordinate system just defined 
(denoted by x e ) can be transformed to coordinates in the standard Apollo 
coordinate system (denoted by Xp) in the following manner: 


X p - (-A Zo )i(-cp d ) 2 (-A 0 ) 3 x e , 

where 

is the longitude of the launch site measured positive west 

cp^ is the geodetic latitude of the launch site 

A z is the launch azimuth of firing direction and is measured 
0 clockwise from a line perpendicular to the plumbline at the 
launch site and pointing northward parallel to the launch 
site meridian around to the firing direction. 

Note that the inverse transformation to obtain equatorial coordinates from 
plumbline coordinates is easily remembered due to its mnemonic form. That 
is. 


x e = <V3(9 d )2(A Zo )iV 
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Also, the transformation matrix, [A] = (A Q ) 3 ($3)2 (A^ ) x , remains constant 
with time. 

For reentry work, a more interesting coordinate system might be a 
coordinate system with the origin at the reentry point with its x-axis 
pointed radially away from the earth- s apparent gravitational center, 
its z-axis parallel to a projection of the instantaneous velocity vector 
in the local horizontal plane, and the y-axis finishing a right-handed 
coordinate system* A transformation from this coordinate system (denoted 
by x r ,y r , z r) to t ^ le equatorial coordinate system can be written as 



where 


A^ is the longitude of the reentry point measured positive west, 

cp_i is the geodetic latitude at the reentry point, 
r 

A^ is the instantaneous reentry flight direction measured as is &z 0 > 
R r is the reentry point altitude. 

Thus, the equations of motion for a point mass in a three-dimensional 
inertial Cartesian coordinate system whose orientation to the standard 
equatorial coordinate system is known can be written as will be shown 
in the following paragraphs. 


2 . Acceleration of a Point Mass in an Earth-Centered Arbitrarily 
Oriented Inertial Coordinate System 

In an arbitrarily oriented earth-centered Cartesian coordinate 
system, the differential equations to be used for simulating the motion 
of a point mass subject to gravitational, aerodynamic, and thrusting 
forces can be written as follows: 


x = (F/m)u + L±Ji 

m 


(GM)x 

R3 


g(x). 


( 1 ) 
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In the above vector differential equation, the first term on the right 
is the thrusting term, where F is the instantaneous magnitude of the 
thrust, m is the instantaneous mass, and u is a unit vector which is 
considered to be a control variable for the direction of the thrust. 

In the second term, L is a vector of aerodynamic forces due to lift 
and 5 is a vector of aerodynamic forces due to drag. The determination 
and control of the thrusting force and the aerodynamic forces will be 
explained in the paragraphs to follow. The third term in equation (1) 
is the gravitational force term for a spherical central body, where GM 
is the gravitational constant and R is the magnitude of the position 
vector of the point mass under consideration. The fourth term which is 
considered to be a function of the position vector is a symbolic repre- 
sentation of corrections to the gravitational force expression which 
are needed when a nonspherical shape for the central body is considered. 
For simplicity in the rest of this development, the fourth term on the 
right-hand side of equation (1) will be ignored, and no further explana- 
tion of the derivation of the third term will be attempted since this 
expression is a very standard representation of the gravitational force 
due to a spherical central body. 

To begin an explanation of the first term (the term due to thrust) 
in equation (1), the thrust force F will be assumed to be constant for 
vacuum flight and to have the following form for atmospheric flight: 


F « F 4- A (P - P) , 
s e o 75 


( 2 ) 


where 

F is the vehicle’s thrust measured at sea level (P = P_) 
s N o' 

A e is the exit area of the engine 

P Q ,P are the atmospheric pressure at sea level and at any 
altitude, respectively. 

For atmospheric flight, the mass flow rate m is given by the relation 


m 



g o I sp 


( 3 ) 


where Isp s is a constant (measured to correspond with F g ) indicating the 
efficiency of a particular set of engines and g D = GM/R§ is a constant 
giving the acceleration due to gravity at the assumed R Q (radius of the 
spherical central body). For vacuum flight the mass flow rate is 
assumed to be the same as for atmospheric flight and the force F v and 
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efficiency I S p v for vacuum flight are related to F g and I S p g as follows: 

F = F + A P 
v s e o 


<VV ^P. 


In either case, the mass at any time t referenced to an initial time t Q 
is given by 


m (t) = m(t D ) + m(t - t Q ) . 


( 4 ) 


To explain how the unit vector u and the aerodynamic force vectors 
L and D are to be determined, a missile-fixed or vehicle-fixed coordinate 
system (denoted by x m ,y m ,z m ) must be introduced. In this coordinate 
system, u will be denoted by and another unit vector a m (a m is the 
vehicle axis which when aligned with the air stream allows no aero- 
dynamic lifting force to be produced) will be introduced. Figure 1 
illustrates the angles to be used to locate and a m in the missile- 
fixed coordinate system. 



Figure 1. Angles used to Locate u m and a m 
in the Missile-Fixed Coordinate System 
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In the missile-fixed coordinate system, the angles cp u and \jr u are meas- 
ured on the ground as a function of trajectory time or burning time to 
locate u m , which is the thrust direction produced by the unswiveled 
engines. Similarly, the angles cp a and are measured on the ground 
as a function of trajectory time to locate a m . These angles locating 
u m and a m in a missile-fixed coordinate system need to be determined 
only once before the initial flight of each vehicle. Then the angle a D 
shown in figure 1 can also be computed as a function of trajectory time. 
The angle a 0 and the unit vector a or a m will be used later in the 
determination of L and D, but for now only the unit vector u or u m is 
of interest. 

As shown in figure 2, the angles Xp and Xy will.be used to locate 
u or u m in the original inertial Cartesian coordinate system (x,y,z). 

The prime coordinate system (x 1 ,y f ,z f ) with the x f coordinate aligned 
with u will be shown again in figure 5 and used to help define the aero- 
dynamic forces. It will be assumed that the vehicle 1 s control system 
acts instantaneously so that the angles Xp and Xy, and thus the direction 
of u in the inertial coordinate system, can be changed instantaneously. 

As the angle Xp and Xy or the direction of u change with respect to the 
inertial coordinate system shown in figure 2, the orientation of the 


X 



Figure 2. Angles Used to Locate u in the Inertial Coordinate System 
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missile-fixed coordinate system will also be changed with respect to 
the inertial coordinate system, since the angles cp u and cannot be 
changed by the vehicle’s control system. 

Often the vehicle’s control system may be given a time history of 
the angles Xp and Xy to follow. In this case, u must be computed in 
terms of X p and Xy so that it can be used in equation (1) for simula- 
tion purposes. As shown in figure 2, 


( 6 ) 


U 1 


u s 

= 

_ u 3_ 



cos X cos X 
y p 

sin X y 

cos X sin X 


If optimization theory is used to obtain u, then the components of u in 
the inertial coordinate system may be determined directly by the optimiza- 
tion theory. In this case, for the information of the vehicle’s control 
system, the angles X p and Xy can be computed in terms of the components 
of u as follows: 

X = arctan (-=—=— i) (- § < X < f) (7) 

y Vuf + u|/ 2 y 2 

Xp = arctan (u!/u 3 ) (0 ^ X^ ^ 2 it). (8) 

This completes the present discussion of u which appears in the 
first term of equation (1). In section IV- C, how u is determined by 
optimization theory is discussed. 

Now to explain how the expressions for L and 5 are obtained in the 
second term of equation (1), the unit vector a or % introduced in fig- 
ure 1 will be used. First, it will be assumed that a ^ u so that the 

unit vectors a and u define a plane in the reference coordinate system 

(see figure 3). 


15 



X 



Figure 3. The Plane Formed by u and a in the Inertial Coordinate System 


In a later section, the simplifying assumption a = u is made and the 
results of this assumption are treated. As mentioned previously, the 
unit vector a is defined to lie along the instantaneous zero lift axis 
of the vehicle. Thus, if a and the relative velocity vector V r of the 
point mass are aligned, L will be zero. The relative velocity vector in 
the reference coordinate system is given by 


V r =x-u>xx-W, (9) 

where u) is the rotational velocity vector of the central body in the 
reference coordinate system and W is an arbitrary vector function describ- 
ing winds in the reference coordinate system. Since a is a unit vector 
in the zero lift direction, a lift force (denoted by L) will be produced 
when there is an angle a between a and V r as shown in figure 4. Figure 4 

shows that the lifting force L which occurs when there is an angle a 

between a and V r lies in the plane formed by the vectors 5 and V r . 
Naturally, the vector D is in a direction opposite to V r . 

To explain further how the unit vector a is to be defined, it will 
be assumed that the directions of u and a in the missile-fixed coordinate 

system are determined as shown in figure 1. It will also be assumed that 

these directions in the missile-fixed coordinate system depend only on 
the vehicle characteristics such as the position of the center of gravity 
and the center of pressure as functions of trajectory time. Thus, the 
angle a 0 between u and a (as shown in figures 1 and 3) can be determined 
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as a function of trajectory time. Since 0C o depends on the vehicle 
characteristics, it cannot be used as a control variable and varied by 
the vehicle 1 s control system. Although the angle Ct Q is not considered 
to be a control variable, it will be assumed that the vehicle 1 s engines 
can roll the vehicle about u (as shown in figure 3) so that the orienta- 
tion of the plane formed by the vectors a and u is a control variable. 
Also, as the direction of u is changed (with the control variables 
Xp and Xy), the direction of a will change because the angle between 
a and u (denoted by a Q ) is predetermined and cannot be varied. 

To write expressions for the components of a in the reference coordin- 
ate system, a new coordinate system whose orientation changes with tra- 
jectory time (denoted by x 1 ,y f ,z f ) will be introduced. This new 
coordinate system is obtained by performing in succession_the Euler angle 
rotations (Xp) 2 and (Xy ) 3 so that the x f axis lies along u as shown in 
figure 2. In the prime coordinate system, the angles a Q and X R (which 
orient the plane formed by u and a) can be depicted as in figure 5. 

The symbol a 1 will be used to denote the unit vector a with its com- 
ponents in the prime coordinate system. These components (see figure 5) 
are given by the following equation: 


a' 


cos a Q 

sin a o cos X R . 

s in Q!q s m X R 


( 10 ) 
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y' 


Figure 5. The Orientation of the Plane Formed by u and a 


Then, a with components in the reference coordinate system is obtained 
as follows: 


a = (X p ) 2 (-X y ) 3 £' . 


(ID 


Notice that the successive Euler angle rotations (-Xy) 3 and (Xp) 2 are 
the inverse of the Euler angle rotations used to obtain the prime 
coordinate system. With a determined by equation (11) in terms of 
Xp, Xy, Xr, and a Q , unit vectors in the direction of L and D can be 
written as: 


(V X a) X v 
v r y r 

I (V r X I) X V r 


u 2 



( 12 ) 


(13) 


Figure 4 will aid the reader in visualization of the construction of 
the vectors u x and u 2 which are written in the reference coordinate 
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system (since a and V r are in the reference coordinate system). Now 
the following expressions for L and D can be written: 


L = - p A (V 
2 H r K x 


V ) C.Ui 
r L ± 



P A (V 
H r v r 


V C D 5 2> 


( 15 ) 


where 

p is the atmospheric density 

A is the reference area of the vehicle 
r 

C is the lift coefficient of the vehicle 

i-i 

is the drag coefficient of the vehicle. 

The atmospheric density p is usually considered to be a function of 
instantaneous altitude and the quantities A r , C^, and C D must be meas- 
ured on the ground. In the measurement of A r , C^, and C-q, a vehicle- 
fixed coordinate system is used (denoted by x v ,y v ,Zv)* For simplicity, 
it will be assumed that the instantaneous orientation of this vehicle- 
fixed coordinate system is defined so that the x v -axis lies along a, 
the z v -axis defines the plane of the vehicle where the relative veloc- 
ity vector V r would produce the maximum values for Cl and Cj) for a given 
angle Ot between and a, and the y v -axis completes a right-handed 
coordinate system. This coordinate system, with u and two possible 
locations for is depicted in figure 6. 



Figure 6. Coordinate System Used to Measure 
the Lift and Drag Coefficients 
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The reference area Ar will be taken to be the area of the(y v ,z v )- 
plane obscured by the vehicle if a sighting along the x v axis is 
made at an infinite distance from the origin. If V r lies in the 
(x v ,z v )-plane as depicted by Vj.^ in figure 6, then the lift coefficient 

Cp and the drag coefficient C D can be measured to_depend only on the 
instantaneous value of the angle a between a and V r and the Mach 
number M of the vehicle. If V r is not in the (x v , z v )-plane as depicted 
by V r2 , then Cp and Cp must also depend on the angle Ctj. as shown in 

figure 6. As mentioned previously, the orientation of u with respect 
to a vehicle-fixed coordinate system Can be determined or measured 
before the vehicle is flown. Thus, the angles CLq and Px can be deter- 
mined using the angles measured in figure 1. The angles Ct 0 and p^. are 
needed in relating the vehicle-fixed coordinate system to the original 
reference coordinate system so that the angles CL and Or can be computed 
and used to obtain corresponding values for Cp and Cp. Figure 6 shows 
that the following equation will determine a for either position of % 
(i.e., V ri or V r2 ). r 


(a X V ) X a 
r 7 


a = arc tan 


(a x V r ) X aj ■ 


(V. • l) 


(-it s a £ it), 


( 16 ) 


To determine 0!r, the angle p 2 shown in figure 6 must be determined. The 
following equation can be used to obtain p 2 : 


P 2 = arc tan 


[a X (u X i) ] . (V r 
(u X a) • (V r X a) 



(0 S p 2 S 2it). 


(17) 


Then is given by 


Or - Pi “ Ps. 


(18) 
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With a and a R computed as shown in equations (16) and (18), values for 
and C]) can be obtained. Thus, enough information is available to 
compute L and D as given in equations (14) and (15). This completes a 
general explanation of the equations of motion given by equation (1). 

To re-emphasize the concepts discussed previously, an outline will 
now be given of the order of the computations to be performed to obtain 
a value for x (as given by equation (1)) at a particular instant of 
trajectory time. At a particular instant of trajectory time (t), the 
position vector (x) of the point mass (m) is known along with the veloc- 
ity vector x. Thus, the altitude (h) of the point mass can be computed 
as follows: 


h = R - R q , (19) 

where R = n/x • x and R 0 is the radius of the spherical central body. 

With a value of h as obtained from equation (19), atmospheric tables can 
be used to obtain corresponding values of pressure (P), density (p), 
and velocity of sound (V s ). Since the gravitational constant of the 
central body (GM) is assumed known, the third term in equation (1) due 
to gravity can be computed immediately. To compute values for the first 
term in equation (1), equation (2) is used to obtain a value for F, since 
A e , F s , and P Q are known constants and P has been obtained as described. 
Also, m and u are assumed to be known at each instant of trajectory time. 
The unit vector u may be known directly in terms of its three components, 
or it may be known in terms of the angles Xp and Xy, in which case equa- 
tion (6) is used to compute u. Then, the first term in equation (1) can 
be computed. To compute values for the second term in equation (1), V r 
must be computed by equation (9). The vector constant w needed in equa- 
tion (9) and the wind function W are assumed known. With a value for 
the vector V r available, the following equation can be used to compute 
Mach number M: 


M = 



( 20 ) 


Also, equations (16) and (18) can be used to compute 0C and a R if values 
of a are available. The unit vector a may be assumed to be known 
directly in terms of its components, or it can be specified in terms of 
Xp, Xy, a Q , and X R as shown in equations (10) and (11). With values 
for a, 0J R , and M obtained, a curve fit or table look-up can be used to 


21 



obtain values of C L and C D , which depend on a , O^, and M. Then equa- 
tions (14) and (15) can be used to compute L and D since the quantities 
needed are now all available. Thus, the third term in equation (1) can 
be computed to complete the computation of equation (1). 

The previous discussion shows how values of x can be obtained at 
each instant of trajectory time when the control variables Xp, Xy, and 
Xj£ are known. Thus, numerical integration of equation (1) could be used 
to obtain a trajectory for prespecified control functions Xp, x , and X R . 

Also, optimization theory could be applied directly to the equations 
of motion (equation (1)) in its present form to determine the correct 
time histories for u and a, but the resulting computer program would be 
much too complicated to consider for guidance purposes. 

For this reason, some simplifying assumptions and approximations to 
these differential equations will now be made that will result in a much 
more compact and efficient computer program for the application of 
optimization theory. 


3. Simplifications to General Three-Dimensional Equations of 
Motion with Atmospheric Effects 

The original equations of motion (equation (1)) are completely 
accurate representations of the ascent phase of the shuttle vehicle’s 
flight. For the vacuum phase of a shuttle vehicle’s flight, an accurate 
representation is obtained by leaving out the term 


L + D 
m 

and for the reentry phase, an accurate representation is obtained by 
leaving out the term 


(F/m) u. 


The only simplifications to be made to equation (1) before optimization 
theory is applied are simplifications of the expressions for L and D 
which will be used in both the ascent and the reentry phases. Also, 
the expression g(x) accounting for oblateness effects will be neglected 
in all three phases, at least in illustrating the application of opti- 
mization theory to the equations of motion (equation (1)). 
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The first simplification to be made in the development of new 
expressions for L and D is to assume that a - CL For most vehicles, 
a Q (the angle between u and a as shown in figure 2) is not very large 
because of the usual aerodynamic symmetry about the thrust vector 
direction. Thus, the assumption a = u is not very restrictive as will 
be seen. With the assumption that a = u, the vehicle- fixed coordinate 
system shown in figure 6 does not have to be introduced. Instead V r 
can be pictured immediately in the prime coordinate system of figure 8. 


x' 



Then, the new expressions for L and D become 


L = t p A (V 
2 r v r 


V ) 
r 


(V X u) X V 
_j£ _ r 

L| (V r X u) X V r H 


( 21 ) 


D = •=■ p A (V 
2 r r 


V ) 

x 



( 22 ) 
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Now the coefficients and C D will depend only on Mach number 
and the angle OC defined by the equation 


(V X u) X V 
r r 

■I (V r XU) X V r | 


a = arctan 1 


• u 


L V r 


(defined for -jt < OC < rt) . (23) 


Also, the following additional approximating equations for and 
Cq will be used. 


C « C, sin a (24) 

L a 

S “ C A + 2 " -cos a) (25) 


where Cqq,, C A , and T) are all now considered to depend only on Mach 
number. These approximations have been found to closely represent most 
C L and Cp tables in the range of Mach and alpha of interest and are 
used mainly because of their simplifcation of the optimization equations. 

Substituting the above approximations into the expressions for L 
and D gives 


L - 2 P A r (V r • V r ) < (C^ sin a) 


(V X u) X V 
r r 

(V r Xu) XV. 


(26) 


= \ P A r (V r • V.) |[C A +. 2rj 0^(1 - cos a)] } 


v r H 


(27) 


24 



Notice that 


| (V r X u) X Vj = |vj 2 sin a - (V. • v_) sin a 

because (X is the angle between V r and u. Also, 

V r • u = ] V r | cos a. 

Thus , 

1 = \ p A r % [( \ X “> X V (28) 

6 = - \ P A r [|V r l (C A + 2n c[ a ) - 2tj (V. • u)] V r . (29) 

Now the following identity will be used: 

(V r x u) X V. = ! V r I s u - (V r . u) V r . 

Then the final form of the lift and drag equations to be used is 

L = \ P A r % C I V r 1 2 u - (V r • a) v r ] (30) 

5 = - | p A r t|v r l (C A + 2 n C^) - 2n (V r • u)] V r . (31) 

These are the expressions for L and D that will be used in equa- 
tion (1) for the application of optimization theory. At this point, it 
will be mentioned again that equation (1) can be integrated numerically 
to yield a trajectory (without any optimization considerations) if a = u 
is known directly in terms of its components or if X p and X y are known 
and equation (6) is used to obtain u. Appendix I is a computer program 
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listing and a sample printout which uses the above equations for L and 
D and assumes the following time history for a - u: 


_ _ 5 o + 3o(t - t o ) 

1% + S o (t " fc o> 


(32) 


* 

The values, u Q , u Q , and the flight time, t^, are parameters in this 
program which can be varied to satisfy specified end conditions and 
minimize desired quantities* In determining how to vary the parameters 
u Q and u Q to satisfy the end conditions, variational differential equa- 
tions for the matrices 


* 

dx(t) c)x(t) 

_ 5 _ 5 

du Q du Q 


Mltl 5 and 


were developed, as explained in reference 14, and an iterative algorithm 
constructed* The program is mentioned at this time so that the interested 
reader can examine the listing to see a sample set of data and how the 
atmospheric and aerodynamic quantities are determined for the numerical 
integration. Other references to this listing and further explanations 
of the technique will also be given in some of the following sections* 


4* Alternate Determination of a and u 

Often it is more desirable to determine a and u in terms of angle 
of attack a, bank angle (3, fa (to be defined), and a Q , mentioned pre- 
viously, instead of using the inertial angles X p , X y , X R , and a Q . How 
this can be done will be explained at this point before discussing the 

application of optimization theory to equation (1). To begin, a new 

coordinate system denoted by x", y", z" will be defined as shown in 
figure 8. Note that, in this coordinate system, the origin is at the 
point mass (m) and r is the position vector of the point mass (i.e., 
r = x in the reference coordinate system). As can be seen from 
figure 8, the V r vector is defined to lie along the z" axis, the 

V r X r vector lies along the y" axis, and the (V r X r) X V r vector lies 

along the x" axis. Thus, the x" coordinate system is transformed to 
the original x,y,z coordinate system as follows: 


x = 


" ( ^r 

/~s 

1 u 

X 

XV r 

V X r 
r 

V -i 
r 

1 <*r 

X r) 

X V 

r 

Iv x 51 

i r 1 

| v J 

1 r 


x 


[B]S n . (33) 
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The inverse transformation is given by 


a" = [B] T a 

and then 


a = arc tan 



(defined for 0 < a < Jt) 


(34) 


S3 = arctan [(a'^/a^)] (defined for 0 s p s 2rt) . (35) 

When ii is not the same as a and is allowed to be oriented with respect 
to a and V r , then another coordinate system must be defined (denoted by 
x'" ,y'" ,z’" ). To visualize this orientation see figure 9. 



Figure 9. The Location of u with Respect to a 
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The transformation from a vector x"' in this coordinate system to the 
original x coordinate system is obtained as follows: 


v r X a 

a X (V r 

X a) -] 

|V r Xa| 

1 ® x (V r 

x S)|J 


[C]x m . 


Thus , 


(36) 


u = [C]u'". 


Therefore, 


a = [c] 


cos a„ 


-sin Q! d sin p r 

-sin a Q cos p r 


Then the inverse transformation is given by 


u"' = [C] T u 


and thus 


(u^’ ) 2 + (u^' ) 2 

a Q = arctan (defined for 0 < a o < it) (37) 

l 


P r = arctan [u£ /u^' ] (defined for 0 I 2x). (38) 
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The preceding equations show how to compute a and u with the 
angles a, p, (X Q9 and P r or conversely how to compute a, p, a Q , and P r 
when the unit vectors a and u are given. If it is assumed that a = u, 
then the triple prime coordinate system does not have to be introduced, 
and only the angles OL and p are needed to compute a and thus u. Con- 
versely, OC and p can be computed when u = a is known or given. 

In the next section, optimization theory will be applied to the 
equations of motion (equation (1)) with L and D given by the simplified 
expressions which resulted from assuming a = u. One of the results of 
the application of optimization theory will be a time history for u, 
and thus a that will produce an optimum trajectory to the desired 
destination. This time history (u = a as a function of time) can be 
used to compute Xp and Xy or a and p as functions of time for guidance 
purposes. 


C. Application of Optimization Theory to Equations of Motion 
1. Equations of Motion 

The simplified equations of motion discussed in the previous sec- 
tion will be used here to illustrate the application of optimization 
theory. These equations with the simplified L and D should adequately 
represent the motion of most space shuttle type of vehicles during 
atmospheric flight. If it is found that they are not adequate, the 
reader should remember that they are used here only to illustrate the 
approach, and whatever corrections are needed can be added to the equa- 
tions without changing the philosophy of the approach. 

The simplified equations for L and D can be substituted into equa- 
tion (1) of the previous section to give the following form for the 
equations of motion: 


GMx 


x = (F/m) u - + g(x)+ (pA r /2m) [V r • V r )u - (V r • u)Vj 


- C(V r • V r ) l/2 (C A + 2r, c£ a ) - 2r\ (V. • u) ] . 


(39) 


In equation (39), F is given by equation (2) and m by equations (3) 
and (4). Typical values of the constants needed to use equations (2), 
(3), and (4) can be found in the computer program listing and sample 
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case printout of Appendix I. The constant GM needed to compute the 
second term of equation (39) and the form of the equations for g(x) 
can also be found in Appendix I, The atmospheric subroutine shown in 
Appendix I shows how the atmospheric density p, the atmospheric pres- 
sure P, and the velocity of sound V s are obtained as functions of alti- 
tude for use in equation (39). With a value for V s , equation (20) can 
be used to compute the Mach number. Then the aerodynamic subroutine of 
Appendix I shows how values of Cl C^, and q are obtained as functions 
of Mach number. If it is noted that V r is given by equation (9) and 
that the instantaneous radius of the oblate earth 1 s surface is deter- 
mined in the differential equation subroutine of Appendix I, then the 
computation of all the terms in equation (39) has been explained. 

A transition to state vector notation will be made to facilitate 
the explanation of the application of optimization theory to equation 
(39). To do this, let 


— - 


— * “ 



X 

x £ 


y 

x 3 


z 

x 4 

= 

X 

x 5 


• 

y 

x 6 


• 

Z 

_ x 7_ 


m 


(40) 


Then 


— — 


— — 


— 

— 



X 


fl( x l, x 2, 

X T ) 

x 2 


y 


f2(Xl, x 2, 

• • ■ 5 ^-7 ) 

x 3 


# 

z 


• 


x 4 

= 

X 

= 

• 


*5 1 


y 


• 


*6 ! 


z 


. 


x 7 


m 


f 7 ( X 1 j x 2 j 

• • • j *^7 ) 


(41) 
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where 


fi = x 4 
f 2 = X 5 
f 3 = X 6 


r f 4i 

f 5 

‘-fe - * 


GM 


= (F/x 7 )u - 


x 2 

I-X3 


(x® + xj + x§) 3 / 2 


+ g(xi,x 2 ,x 3 ) 


P A v 


Ct^ [(V • V)u - (V • u)V ] 


2x^\ u La L ' r r 


t(V r • V r ) 1 / 2 (C A + 2r) cf a ) - 2 t] c£ a (V r . u)]V 


(42) 


f 7 = m (a constant developed in equation (3)). 

Note that, in the previous expressions for f 4 , f 5 , and f 6 ; F 
depends on x 1 ,x 2 ,x 3 ; p depends on x 1 ,x 2 ,x 3 ; V r depends on xi,x 2 ,x 3 ,x 4 , 
x 5 ,k 6 ; and the aerodynamic coefficients C^q,, Ca, and r\ all depend on 
x i> x 2 , x 3 > x 4 , x 55 x 6 " This dependence can also be seen by referring to 
the differential equation subroutine of Appendix I. The preceding 
transition to state vector notation now allows the equations of motion 
to be written in standard first order form which is usually seen when 
optimization theory is developed. That is, 


X = f (X, u) , 


(43) 


where X, f, and X are all vectors of dimension 7, and u is a vector of 
dimension 3. 

Before applying optimization theory to this system of equations, 
the optimization criteria and the boundary conditions must be discussed. 
Such a discussion appears in the next two sections. 
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2. Definition of Optimization Criteria 

The optimization criteria selected for discussion in this section 
were chosen mainly for illustration, and some approximations were 
made to keep the equations as simple as possible while maintaining as 
much physical realism as possible. Again it must be emphasized that, 
if this particular selection of an optimization criteria is found to 
be unsatisfactory, then whatever criteria that can be selected which 
are satisfactory can be used without changing the basic philosophy of 
the MASCOT technique. 

The selected criteria are given as: 


t 


f 


Xs(t f ) - kix 7 (t f ) 4- k 2 



Q dt 


+ k^ 


f D 2 dt . 
t 

o 


(44) 


That is, x 8 (tf) is to be minimized where the constants or weighting 
factors ki, k 2 , and k 3 can be selected to achieve different ratios of 
values among the three terms in the expression for x 8 (tf). For example, 
in reentry, could be zero with k 2 and k 3 equal to .5, or in ascent, 
kx could be minus one and k 2 and k 3 zero. In equation (44) Q, the con- 
vective heating rate per unit area of the vehicle, is given by 


Q = 


\To 


W 2 



(45) 


where e is a constant and a is the radius of curvature of the vehicle 1 s 
nose. Values for these constants are given in Appendix I. Also, in 
equation (44) 


D* 


- 2 - 2 

L + D 


m 


The following steps 
From the simplified 
can be seen that 


are needed to derive the exact expression for D . 
equation for L and D (equations (30) and (31)), it 


L| 2 = L * L = (\ p A r C La ) 2 [(V r • V r ) 2 - (V r • V r )(V r • u) 2 ] (46) 
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|Dl 2 =D.D= (|pA r ) 2 [(V r • V 1/2 (C a + 2r] C^) 

- 2f| C 2 (V • u)] 2 (V • V ) 

a v r 

- <? P V 2[ <\ • V< C A + 2 1 <%/ 

' 4 1 c f a • V 1 ' 2 ^ • " )<C A + 24 ^a> 

+ 4 0l ^Q Z <V r • u) 2 ](V r • V r ). (47) 

Now note that (V r * u) = |v r |cos OL so that 

(V r • u) 2 = (V r • v r ) cos 2 a = (V r • V r )(l - s in 2 CC) . (48) 

In equation (48), sin 2 a will be approximated by 2(1 - cos Q!) to give 

(V. • u ) 2 = (V. • V r )[l - 2(1 - cos O')] 

= (V r • V r )[2 cos a - 1] (49) 

=2 (V r • V r ) 1 / 2 (V r • u) - (V r • V r ). 

The substitution of the final form of equation (49) into equations (46) 
and (47) gives 

|£| 2 « (j p A r ) 2 (2C 2 a )[(V r • V r ) 2 - (V r • V r ) 3/2 (V r • u)] (50) 

|D| 2 - (| P A r> 2 t^ r ' V r )(C| + 4n C A C^) 

_ 4T 1 C A c L a (V r • u)(V r • V r ) l/2 HV r • V r ). (5D 
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Adding equations (50) and (51) and dividing by m 2 gives the final form 
to be used for D 2 a That is. 


D 2 - (pA r /2m) 2 [(V r • V r ) 2 (C 2 + At) C a + 2C^) 

- 2(V r • 3)(V r • V r ) 3/2 (C^ + 2i) C a 0^)]. (52) 

With the preceding explanation of Q and D 2 and using equation (44), 
it can be seen that 


x 8 = k x x T + k 2 Q + k 3 D 2 = f 8 (x 1 ,x 2 , . . . ,x 7 ) . 


(53) 


Thus, x 8 (t) as defined by equation (53) is a new variable which is 
added to the seven variables of the preceding section so that x 8 (tf) 
is the quantity selected to be minimized at the final time with 
x 8 (t 0 ) = ki m 0 * 

3. Definition and Explanation of Some Example Boundary or 
Mission Conditions 

The previous two sections have put the illustrative problem being 
solved in this report in the form: 

(a) X(t ) and t Q are given 

(b) X(t) is defined by X(t) = f[X(t),u(t)] (54) 

where t Q ^ t ^ tf 

(c) x 8 (t£) is to be minimized by the choice of Q(t) * 

In equation (54), X is now an eight-dimens ional vector and u is a three- 
dimensional vector. The only thing missing from equation (54) is some 
additional constraints which are usually placed on some of the values 
of x x (tf), x 2 (tf), x 7 (tf). These functional constraints are called 

the boundary (or end or mission) conditions. 
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The first set of boundary conditions to be discussed is an example 
set of reentry boundary conditions. For example, on a reentry trajec- 
tory a desired value for the magnitude of the radius vector, Mach 
number, the path angle, down range, and cross range could be specified 
for the end of the trajectory. In functional form the preceding speci- 
fications result in five equations as follows: 


Q — 


2 2 , 2 
x x 4* x 2 + x 3 


-1 = 0 


G„ - 


V r 


V 2 M 2 
s c 


-1 = 0 


(54) 


(55) 


** / 2 2 ~ ~ 
R c VX 4 + X 5 + X 6 


- COS ■a =0 
C 


(56) 


G4 “ 


r *1 1 

x 2 I • 

x 3 J 


X 


X 


LI o| 


tan %R 


x d x ; ) -1 


x X (x X x ) I -I 
O V o o' 1 



(57) 


G5 



x i ' 
x 2 

x 3 J 



(58) 


In all of these equations, the values of xi, x 2 , X3, X4, X5, and x& are 
used at the final time. In equation (54), R c is the desired value of 
the magnitude of the radius vector at the final time, and when G x =0, 
it can be seen that x 2 + x§ + will be equal to R 2 . Similarly, in 

equation (55), M c is the desired value of Mach number at the final time; 
in equation (56), 6 C is the desired value of the angle between the 
inertial position vector and the inertial velocity vector (path angle); 
in equation (57), cpj)R is the desired value of the down range angle; 
and in equation (58), cp CR is the desired value of the cross range angle. 
In order to better understand equations (57) and (58), figure 10 is used. 
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Figure 10. A Definition of the Down Range and Cross Range Angles 

In this figure a coordinate system is defined using the position and 
velocity vector at the initial time (x Q and x 0 ), and then the angles 
cpDR and cp^ are used to locate the final position vector x^ with 
respect to the coordinate system defined by x Q and x Q . The figure 
shows that 


x X (x X x ) -] 

o o o 

L l*„ x (So X S„)|- 


tan 'Pdr ■ 


O' -J 


8in ‘PoR 


f x o X x o 


<- | x Q X x Q 1 -J 


(59) 


(60) 


Since |xf| is to be equal to Rc, it can be seen that equations (57) and 
(58) follow from equations (59) and (60). 
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The computer program for computing equations (54), (55), (56), (57), 
and (58) is called subroutine BOUND in Appendix I # Notice that all of 
the equations are normalized so that G| + G§ + G§ + Gf 4* G§ will have 
a magnitude near the value of one* The preceding sum of squares is 
usually used as an error check to determine when the boundary conditions 
are satisfied. That is, it must be less than some tolerance (for 
example, .5 x 1CT 10 ). 

With the preceding definition of an example set of boundary condi- 
tions, the complete trajectory optimization problem being solved in this 
paper can now be stated in functional notation* That is, 


(a) 

X(t Q ) and t Q are given 



(b) 

X(t) is defined by X(t) =* 
where t Q ^ t ^ t^ 

f[X(t),u(t)] 

(61) 

(c) 

Xq ( t f ) is to be minimized 
made equal to zero by the 

and G[Xf(tf),tf] is to be 
choice of u(t). 



Appendix I is a computer program listing which solves the problem stated 
in equation (61) with the assumption that 


5 o + - fc o> 

“o + "oCt " t o'> 

where the vector constants u Q and u Q are assumed to be unknown param- 
eters which the program must select along with tf to satisfy condition 
(c) of equation (61). This approach to trajectory optimization is 
called parameter optimization particularly when the time interval t Q 
to tf is broken into n segments (at the times t l 5 t 2 , . t n ) and u(t) 
is defined by 


u + u (t - t ) 
o o o 

u(t) = ; — for t 0 S t S ti 

|u 0 + S o(t - to) I 

u(ti) + u x (t - t x ) 

u(t) = 7 for t x S t S t 2 

u(ti) + U!(t - t x ) 


u(t n ) + u (t - t n ) 

u(t) = for t n I t s tj 

5 (t n ) + « „(t " t n ) | 


( 62 ) 
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Then the parameters the program must select to satisfy condition (c) of 
equation (61) are u Q , & 0 , u x , u 2j u n , and t f . As stated previously, 

the computer program listed in Appendix I and the sample case are a 
solution of the problem stated by equation (61) using only the param- 
eters u Q , Q 0 , and t£ with the boundary conditions given by equations 
(54), (55), (56), (57), and (58). Before explaining in more detail how 
this computer program works, a set of example boundary conditions for a 
rendezvous mission and an ascent mission will be given, and it will be 
noted here that the computer program of Appendix I has been used to solve 
these problems when the subroutines BOUND and PBOUND are replaced with the 
appropriate new subroutines for ascent and rendezvous. 

For an ascent mission, the magnitude of the radius vector, the magni- 
tude of the velocity vector, the angle between the radius vector and the 
velocity vector, and the angle between the plane containing the position 
and velocity vector and the equatorial plane could be specified. This 
yields the following set of boundary conditions (note that all values of 
xi, x 2 , x 3 , x 4 , x 5 , and xg are used at tf): 


Gn - 


2 

*2 


, 2 

+ X3 


-1 = 0 


K 


(63) 


G s = 


2 2 2 

X 4 + X 5 + X 6 



1 


0 


(64) 


x l x 4 + X 2 X 5 + X 3 X 6 

G 3 cos fl c = 0 


(65) 


G 4 = 


r 

r Xi - 


r x 4 " 

I 

i 

x 2 

X 1 

x 5 

f * 


L X 3 J 


L X 6 J 

J 


r v„ 

c c 


- sin cp„ cos i c = 0. 


( 66 ) 


In equation (63), R c is the specified radius vector magnitude; in equa- 
tion (64), V c is the specified velocity vector magnitude; in equation 
(65), *Q C is the specified angle between the radius vector and the veloc- 
ity vector; and in equation (66), i c is the specified inclination of the 
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flight plane to the equatorial plane. Also, in equation (66), A is a 
unit vector perpendicular to the equatorial plane. Figure 11 shows A 
and a plumbline coordinate system oriented, with respect to the equa- 
torial coordinate system, by the azimuth Az and the geodetic latitude 
cp d of some point such as the launch site. 



y 

Figure 11. The Orientation of the Plumbline Coordinate System 
and the Unit North Vector A 


Figure 11 shows that 


A = 


sin cp d 

cos cp^ sin A z . 

cos cp d cos A z 


(67) 


This completes the explanation of an example set of ascent boundary 
conditions . 
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For a rendezvous mission, let the coordinates of the target vehicle 
(in the same coordinate system as the pursuit vehicle) be given by the 
six-dimensional vector Y, where 


X 


yi 

y 


yz 

2 



X 


y4 

y 


ys 

# 

z 


_ys_ 


Then the boundary conditions 


Gi - 


Gs - 


fi 

Yi 

X2 

y 2 


- 1 


- 1 


G 3 = — 
Y3 


G 4 — 


y 4 


- 1 


g 5 = — 

ys 


G6 


XS 

ys 


can be written in the form: 
« 0 
= 0 
= 0 
= 0 
= 0 
= 0. 


( 68 ) 


(69) 

(70) 

(71) 

(72) 

(73) 

(74) 


Note that in the above equations the values of x x , x 2 , x 3 , x 4 , x 5 , and 
x 6 are all taken at the final time and the values of y x ,y 2 , . ye 
depend only on the final time so that equations (69), (70), . . . , (74) 
form a set of boundary conditions in the form G[X(tf),tf] = 0 .where €. 
is a six-dimensional vector. 

The preceding discussion completes the explanation of all three 
sets of boundary conditions. Now the operation of the computer program 
in Appendix I will be explained in more detail. The starting point will 
be the problem statement of equation (61) with the five -dimens ional 
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G[X(tf),t f ] vector given by equations (54), (55), (56), (57), and (58). 
As stated before, u(t) is assumed to have the form 


u(t) 


u + U (t - t ) 
o o o 

u + U_(t - t ) 
o o v o' 


and the parameters u 0 and u Q (appearing in the u(t) expression) and t£ 
are used to satisfy condition (c) of equation (61). With the preceding 
form for u(t) and the fixed initial conditions of condition (a) in equa- 
tion (61), it can be stated that X(tf) depends only on OC , P, u Q , and t£, 
where OC and (3 are shown in figure 8, and define a unit vector for u Q as 
shown in the explanation following equation (33). That is. 


u = (sin OC cos (3) 


r (V r X x ) X V r 
r o o r o 

-| (V r X X) X v r I 


+ sin a sin p 


r Vf o x x o 
HV ro X Sj-i 


+ cos a 


v r 


v r 


(75) 


It might also seem that X(t^) should depend on the magnitude of u Q , but 
the definition of u(t) as a unit # vector tends to eliminate this possi- 
bility because the magnitude of u Q must also be considered. That is, 
any arbitrary initial magnitude of u Q will produce the same X(tf) that 
an initial magnitude of u Q equal to # one would produce, when OC and |3 are 
kept the same and the magnitude of u Q is adjusted with the ratio of the 
two different u Q magnitudes. 

Now it can be stated that G[X(t£),tf] and x 8 (t£) depend on OC, |3, 
u Q , and t £ . Thus, two new vectors (F and z) can now be defined. That 

is. 


F(z) = F( 0 !,P,u o ,t f ) 


G[X(t f ),t f J 

Wix 8 (t f ) 


(76) 
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where 


a 


z = 


P 

A 

u 

o 


In equation (76), notice that F is a six-dimensional vector and z is 
a six-dimensional vector. The fact that F and z are of the same dimen- 
sion is somewhat coincidental, since, for the following development, 
it is required only that the dimension of z be greater than or equal 
to the dimension of F. Note also the arbitrary weighting factor VS 1 
whose use will be explained later. 


The preceding discussion allows condition (c) of equation (61), 
to be stated in the following form: Determine the parameters a, |3, u, 

and tf to minimize F^ F with a choice of the weighting factor Wi 
which allows G to be zero. Note that F^ means F transpose so that 
F^ F is a scalar, and similarly for G^ G. To determine the vector of 
unknown parameters z, which minimizes F T F, a Taylor series expansion 
of F about an arbitrarily guessed vector z* is made. That is, 


F(z) = F (z*) + (Sf/5z)^(z - z*) + ... . 


(77) 


The notation (dF/dz)* means the matrix of partial derivatives of the 
vector F with respect to the vector z evaluated at the star point. 
Its determination will be explained later. If (Sf/Sz)* in equation 
(77) is not zero, then a vector z can be determined such that 
F(z) - KF(z*) where K is a scalar such that 0 ^ K < 1. Also, if all 
the terms of order greater than one are neglected in equation (77), 
then the difference in this linear approximation to F(z) and KF(z*) 
will depend on z. That is, 


e(z) = F (z) - KF(z*) 

= (1 - K)F(z*) + (dF/dz)*(z - z*). 


(78) 
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T 

Equation (77) shows that e (z)e(z) is a minimum when the vector z 
satisfies the following equation: 


(SF/Sz)^ (dF/dz)* (z - z*) = -(1 - K) 2 (SF/Sz)^ F(z' v ). (79) 

Notice that, if the function F(z) satisfies appropriate continuity condi 
tions, then a value of K close enough to one can be chosen so that the 
value of z obtained from equation (79) will produce a value of F(z), 
using equation (77), that will be equal to KF(z*) within any desired 
degree of accuracy. The computer program of Appendix I uses equation 
(79) and some logic for choosing a sequence of Wi f s and a sequence of 
K* s to determine a vector z which minimizes F(z). Since the program 
listing is self-explanatory, an explanation of the logic is not 
included here. The only other aspect of the listing in Appendix I to 
be explained is the determination of (c)F/c)z*) which appears in equa- 
tions (77), (78), and (79). 

Since the vector F depends on X(tf) and tf in its original form, 
(c)F/dz)* is obtained using the chain rule. That is. 


F[X(t f ),t f ] F[X(t f ),t f ] 

(dF/dz) = ; = — 

Sc£,j3,u Q ,t£ X(t f ) 

(80) 

On the right-hand side of equation (80), it can be seen that 


£x(t,)\ SF[X(t f ),t f ] /St A 

ST-) + 5^ (j5T> 


SF[X(t f ),t f ] 

Sx(t f ) 


and 


SF[X(t f ),t f ] 

bt,. 


are obtained by explicit partial differentiation of the expression for 
F[X(t^),t ; p] with respect to the variables X(tf) and tf. Also, to 
further explain the other two terms on the far right-hand side of 
equation (80), the following definitions are repeated: 
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dt f 



3 a 


0 


dt f 


0 

dtf 

ST 


0 

sr = 



0 




0 

x'O .1; 



1 

: . i i , 




and 

dx(t f ) pax(t f ) ax<t f ) ax<t f ) 

ciz = da ’ bp ’ ’ X( V_ • 


( 81 ) 


(82) 


In equation (82), the chain rule is again applied to the first two 
columns of the matrix for dX(tf)/dz to give 


c>X(t f ) 

SETT. 


ax(t f ) ax(t ) . 

< aa 0 /Sa) + ' - (du 0 /da) 

0 Ou 0 


and 




Sx(t ) . 


lp“ = .+ —3— (Su o /bp) 

0 0U„ 


(83) 


(84) 


Since u 0 , 0!, and (3 are assumed to be independent parameters, du 0 /da and 
du 0 /dp are zero. Equation (75) gives the relation between a, (3, and u 0 
so that 


du 

o 

da 


8 4 ?r>;i v >• yo.'Jtf* 


= COS a *SCOS P 


{• 


(V ro X 5 0 ) X V ro 
L l(V r „ XL) XV, |J 


- sin a 


-|»r 0 N 


+ sin p 


V r X x 
o o 

Iv r . x L|J 


(85) 
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and 


5u 

o 


sin a 1 cos P 


X x 


V r X X H 

1 r n O 1 


- sin p 


V r X x X V r 
L o o r o 

L l(Vr n X 5 0 ) x V, 


( 86 ) 


Now the only undetermined matrices needed for ciF/Pz, as given by equa- 
tion (80), are the matrices dx(tf)/Bu 0 and PX(tf)/c)Q 0 which appear in 
equations (82), (83), and (84). These matrices must be determined by 
numerical integration of the variational equations. Equation (61) 
shows that X(t Q ) and t Q are fixed and X = f(X,u). Thus, 


3x(t o ) 


0 


Bx(t ) 

o' 


= 0 


(87) 


( 88 ) 


and 


d__ 

-ax(t) - 


3f (X,u) 


af cx,u)~ 

(-ax(t) - 


"afcx,G) 

ras(o ' 

dt 

Lau 0 ,a 0 J 


^ 5 o’“o j 


ax 

^“o» V 

T* 

Pu 



The initial condition matrices given by equations (87) and (88) allow 
equation (89) to be integrated numerically from t Q to tf to determine 
the matrices Px(tf)/Pu Q and Px(t^)/Pu Q . Notice that the matrices 


df(X,u) 


df (X,u) 

- aa J ’ 


and 


au 


L ^n>u 


o’ o 


are obtained by explicit partial differentiation of the expressions 
for f(X,u) and u(t). 
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The preceding discussion completes an explanation and brief 
derivation of the basic ideas included in the computer program of 
Appendix I. For the actual expressions of all the partial derivatives 
and any other equations of interest, the reader is referred to Appendix I. 


4. Derivation of Optimization Conditions 

In the previous section of this report, the optimization problem 
described by equation (61) was solved using the parameter tf and the 
parameters u Q and u c appearing in a particular form of u(t). In this 
section, the optimization conditions needed to determine u(t) satisfy- 
ing equation (61) will be derived without assuming a specified form 
for u(t) (i.e., the optimization conditions will determine the form 
for u(t)). The derivation used here is rather informal, and is 
included mainly to give the reader an idea about the basis for opti- 
mization conditions and their form. A more precise derivation for 
similar problems can be found in references 15 and 16. 

The first step in this derivation is the definition of an optimiza- 
tion criteria (J) with the constraints adjoined to J with vectors of 
"Lagrangian" multipliers (p and A). That is, 


J = x 8 (t f ) 


t 


f 


+ p L G[X(t f ),t f ] + 


t 

o 



X(t)-f[X(t),u(t)] 


dt . 


(90) 


In equation (90) p is the vector of multipliers associated with the con- 
straints of condition (c) in equation (61), and A(t) is the vector of 
multipliers associated with the constraints of condition (b) in equa- 
tion (61). Since p is defined to have the same dimension as G and A 
the same dimension as X, it can be seen that J is a scalar. For a 
given set of X(t Q ) and t Q , a minimization of J, a satisfaction of 
conditions (b) of equation (61), and a satisfaction of the condition 
g[X(tf),tf] = 0 by the selection of u(t) will produce a minimum of 
Xa(tf) and thus solve the problem stated in equation (61). 

To minimize J as given in equation (90), with the proper selection 
of a u(t), consider the class C of all u(t) (where t 0 ^ t ^ tf) which 
satisfy the differential equations X = f(X,u) and the boundary condi- 
tions G[X(tf),tf] = 0. Let CT(t) denote a particular control function 
in this class and 6u(t) a small variation in u*(t) such that u(t) equal to 
u*(t) + Su(t) is in C. Then a*(t) will be at least a local minimum 
of J if 5J = 0 and 6 2 J > 0. The computer program in Appendix I can be 
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used to show that the class C is not empty l and that: variations in u* ( t) 
(produced by a 5 u(t)) are also in "C at t f £or < : shmeoofc 4 he r 'u^(i) in r € # ^ ; 

Thus, the computer program >dn Appendixn I ! can be-uSed^to show"' the exist- ?qn * 
ence of a nonempty subclass^ ^ eontaihed. in' frsuchvthat^'every ^element r 5 • 

u(t) in C* is continuous with respect to t and has variations SQ(t) 
which produce continuous variations in the trajectory given by X(t) a 
This is the class of control functions * to? be cons id efaed:! h^ree inr the * ! ! 

examination of &J and S 2 J. Thus, for a particular tl*(t) in C*, a 
small variation &u(t')^ wili gi^e r :r: * ? I - :;o suol vb •• i c r '"' csi 
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c'X(t) 
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In equation ( 91 ) the term 
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•t « X 
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can be integrated by parts/to ^ive ^ 


>i> ? 


j ^ I ?* ^ 


f*‘‘) CiO i 


f ;1 e r? * t ; 


. t r* ; ; -j j cjc: 00.1 oT- £ >f t - • : ? «< :.;>0 I 

i n f T. ■■ ■ i “t p.:- j.- « r : .< rn k v 0 :* bn b 1 . ,: 'o. > ov.r 

:> n r> f 1 or l * ? b ' f;r - .. ( oop.v lo .< d) : v\< .? bno 

r A T 8X = [A t 5X]/. 7 ' f : a t ( t) ; bx(t) . . 1 r ' : ' ; o , -* ,f 5 , v,i 0 * ; m) i 1 ' 

t t 

O ■ . . o ,, .... 

: .- . ( *'• ! ^ ?* x ♦.» ' nr ir)7 - o ?*S .• •' Q J 

•’ . . < ■: ;■> j’ I : . • ■■ • " P; b 1 ' j ?d i '^'>h ! " ).{.'■ ■; j. 

•• ; J ■■ :••• OBJ h;;?B ( n B % '' “ '! ^<i>; ! A ^.upa j c. ’ irj'Vf Op f!) J o’! A t j ; 

Substitution of equation'. (-92|) ; into equationf (91) and nosting | £h§ty sifiOe - c.t 
X(t Q ) is fixed, 5X(t 0 ) is zero gives r'B-S' '• ! r.«ic» s (l)nc' one ??; ■ s?r>; r»’ 
• •• •; S .; .-. . j * ;i ar?i : p • ’ > > } ’ 3 :fO»fT ... ^ t;.* ;;.j; ' ]' 1 :.1 ‘ '• 


' U'j ’■'> 


f-O:, 0 


48 



■I , ' 'I! ‘.f. ! 11 


5J = 4 10,0,0,0,0,0,0,1] + p 


■ T' 


dG[X(t f ),t f J: 

5x(tp 


i* 


+ a (t f > y sx(t f > 


: i -ii ‘ J i ■■■ h'~ • ■ -V » ; O j. V " / 
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+ A T (t) 


of [X(t),u(t )] 

dx(t) 


SX(t) 


t ‘ ; , i :* 


+ A T (t) 


V(t) J* 5 U(tj j ‘ IC - 


(93) 


n o 


i i ' ex 1 i:. % n ,:j j. ! ?w v i\> J :• ^ f * jc 

•j v O '1 * ' " * „ * jH; f j j"> ‘ . r* ; 

From equation (93), the list of necessary’ conditions are obtained by 
requiring J;^.t 5J = 0 f This gives u . . ■ _ ; r . 


■ } i i O : J r i ' O :J.t v/ U -O 


(a) A T (t.) = -[0,0,0,0,0,0,0,1] - p T 


r dG[X(t f ),t f ] 
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(d) A (t) 


~ cif[X(t),5(t)] 

Su(t) 


= 0. 


(94) 


*. :« i b .:> (> i > ) i : nl xux •(”> : • ,! *„ : i.r- 'l xx , ■ . 

Thus, any candidate trajectory for- the .minimization. of J must satisfy 
conditions (a), (b), (c), and (d) of equation (94) # Conditions (a) and 
(b) are called transversality conditions, condition (c) is the defining 
differential equation for A(t), and condition (d) is used to define 
u(t). The next section of this report will explain more about how 
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equation (94) is used to obtain candidate trajectories for the minimiza- 
tion of J. 

At this point S 2 J could be developed and examined to obtain the 
conditions which assure that > 0, but this has already been done 

and explained very well in chapter 6 of reference 16. Thus, it will 
not be done again here, but the reader is urged to examine this material 
to understand the following sufficiency conditions which will be given. 
That is, S 2 J > 0, if 


(a) 


a 2 {A T ffx(t),u(t)i3 

Su(t) du(t) 


is negative definite , 


(95) 


(b) the normality condition is satisfied on a trajectory which 
satisfies equations (94) and (a) above, 

(c) no conjugate points occur on a trajectory which satisfies 
equations (94) and (a) above. 

In equation (95), condition (a) is usually called the strengthened 
Legendre-Clebsch condition and can be combined with condition (d) of 
equation (94) to say that 

H(t) = A T f [X(t) , u(t) ] (96) 


must be a local maximum with respect to u(t) for t Q S't ^ tf. An explana 
tion of the meaning of conditions (b) and (c) of equation (95) will be 
deferred until section 6 on the solution of the boundary value problem 
resulting from an application of the other optimization conditions just 
developed. This completes the derivation and explanation of the optimiza 
tion conditions which will be used in the next section. 


5. Application of Optimization Equations 

The complete set of differential equations in first order form to 
be used for this illustration is given by 
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Xi = x 4 
x 2 = x 5 

&3 = X 6 


X 4 

x 5 

*6 


(F/x 7 )u 


GM 


- x x - 
x 2 

L x 3 J 


^ x i + x l + x §) 


3 / 2 


+ g(xi,x 2 ,x 3 ) 


+ (pA r /2x 7 ) -|Cj^ [ (V r • v r )u - (v r * u) v r ] 

- [(V r • V 1/2 < c a + 2r, C^) - 2,, C^(V r • 5)] V r | 

• « 

X7 = m 


(97) 


x 8 = kim + k 2 (e/sTa)p 1 / 2 | V r | 3 

+ k 3 (pA r /2x 7 ) e [(V r • V r ) 2 (C| + 4n C A C + 2C^) 

- 2(V r .5) (V r • V r ) 3 / 2 (C 2 a + 2r] C a C^)]. 

With these differential equations the Hamiltonian H(t) defined by equa- 
tion (96) can be written as 


GM 


H(t) - ^1X4 + A 2 x 5 + A 3 x 6 4 - [ A4 , X 5 , 7 \ q ] 


(F/x 7 )u 


x i 

x 2 

x 3 


(x 2 + x| + x|) 3 / 2 


+ i(xi,*s, k 3 ) + (pA./2x 7 ) jcjJ (v r • V r )S - (V r - u)VJ (98) 


[<v r • v r )V= ( c A + 2n c£ a ) - 2r) c^(v r . 5)] v r 


u a: 


(equation (98) continued on next page) 
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+ AylTl 


+ A 8 |k x m + k 2 (e//a)p l/2 |v r | 3 + k 3 (pA ; /2x 7 )P [ (V; • V.) 2 (Cj + 4tj C a C^+20 f^) 


-2(V r • S)(V r • V r )3/ 2 (C 2 a + 2r| C a Cf^) ]J- . 


i* •. i i 


i ..y 


To maximize H(t) with respect to u(t), only the terms involving u(t) 
need to be considered. That is, 


H* (t) 


[ Ag , A 6 ] 


(F/x 7 )u + 


(pA r /2x 7 ) 


{' 


V< V r * V u 


(V r • u)V r ] 


+ 2r\ C 


a 



J 



(tey- 






-A 8 k 3 (pA r /2x/[2(V r • S)(V r r V r ) 3 / 2 (Cr + 2 tj C a C* )]. 

0£ 


y) (: 


Then H(t) can be written as 



H(t) = z(t) + H' (t) = z(t) + h T (t) u(t), 


(100) 



A -i ••• . \ i _ X. . A ( .> h 
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where 


z(t) - A 3X4 + A 2 X 5 + A3X5 + [^ 45 X 5 ^ 5 ] 


r Xi 1 

/ GM ; x 2 
* L X3 - 

(x 2 + x| + x 2 ) 3 / 2 


+ g(Xl,X 2 ,X 3 ) 


XT W 2 


-(pA r /2x 7 )(C A + 2r\ Cj^) (V r * V,.) 


+A 7 m + A a [k-jm + k 2 (e/\T if) p 1 / 2 1 V | 3 


+ k 3 (pA r /2x 7 ) 2 (V r • V r ) 2 (C| + 4q C a + 2C 2 a ] 


(101) 


h = [A 4 ,A 5 ,A 6 ] j|(F/x 7 ) + (pA r /2 x 7 ) C La (V r * V r )j p 1 0 

a - ' ! ) 0 0 1 


10 0 




■ * i>. h '•* • i i -i t v« 


+ (pA r /2x 7 )[(2r l C La -l) C La (V r V*) ] 


(102) 


-A 8 k 3 (pA r /2x 7 ) 2(V r • V r ) 3/2 (CL + 2q C A C^) V*. 


,2 s -T 


Sitice u(t) pas been' assumed . to be a unit vector throughout; this report, 
it cari'be wc£€ten as ' " " ’ ’ v ” J v* » „• j 


|p(t)i • 


(103) 


so that p(t) becomes the actual control variable which is used to maxi<- 
mize H(t). -Applying condition '(d) of equation (94) to equation (100) 
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with 


u = - 2 &- 
|p(t)| 


gives 


Bp(t> 


= r t 




(104) 


As can be seen in equation (104), 


dH(t) = 
dp(t) 


0 , 


when p(t) = ±ii(t). Also, condition (a) of equation (95) applied to equa- 
tion (100) will show that 


d g H(t) 
dp(t) Sp(t) 

will be a negative definite matrix only if p(t) = h(t). Thus, the 
optimal control u(t) for t Q ^ t ^ is determined by 

u(t) = 7^7 , (105) 

|h(t)| 

where h(t) is given by equation (102)* As can be seen in equation (102), 
the values of some of the multipliers . ..,A S are needed at each 

t to determine u(t). Thus, the differential equations for these multi- 
pliers given by condition (c) of equation (94) must be used. This 
equation can be written in the form 


A = 



(106) 


where H is given in equation (98). To take the partial derivatives 
indicated in equation (106) and write explicit equations for the A f s 
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is a rather lengthy task and thus will not be included here. The 
Fortran expressions for these differential equations can be found in 
the differential equation subroutine (DEQ) of the computer program 
listing in Appendix II. For a given set of initial conditions (X(t Q ), 
A(t Q ), and t Q ) , this computer program integrates numerically the dif- 
ferential equations given by equations (97) and (106) and uses equation 
(105) to yield at each t (where t Q ^ t ^ tf) a value for X(t), A(t), 
and u(t). The numerical integration of the differential equations stops 
at a given time (tf). At this time tf conditions (a) and (b) of equa- 
tion (94) can be applied to any of the three sets of boundary conditions 
(G’s) described in section (3) to give the transversality conditions 
which complete the boundary value problem. As an example, equations 
(54), (55), (56), (57), and (58) are used, and the result is the equa- 
tions listed in subroutine BOUND of Appendix II. As can be seen in this 
subroutine, or in conditions (a) and (b) of equation (94), values for 
the end condition multipliers (p f s) must be given to evaluate the 
transversality conditions. 

This completes the discussion of the application of the optimiza- 
tion conditions developed in the previous section. To reiterate the 
steps in the calculation of a candidate optimal trajectory, the following 
steps are listed and called equation (107): 


(a) 

(b) 

(c) 


Values of X(t Q ) and t Q are given as fixed values. 

Values of X(t Q ), p, and tf must be known or determined. 
The differential equations 


X(t) f [X(t) , u(t) ] 

MO = 


cif[X(t),Q(t)] 

c)X(t) 


Mt) 


(107) 


with 


u(t) 


h(t) 
|h(t) | 


(where h(t) is defined by equation (102)) must be integrated 
numerically to tf with the initial conditions given in condi- 
tions (a) and (b) . 


(equation (107) continued on next page) 
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(d) The known or determined values, of X(t 0 ), p, and as,jgivep , . t 
- in (b) must cause the*, following, boundary condition to be . 


-.satisfied: : 

G[x(t f >,t f ] == o 


- $ i < 






- -■ .i.i J O ■ ' t 't j A ■ ■ t?/I J 

h ^ rbQ[X(t f ) ,t Jt ? j - ^ t-ac 

i • A A (t f ) » -i[.o j, o , o , ,o , o , o , o.-, i ] ^ p 1 — dxxi ^ 

< -A ■■;■.■ 1 :j . : V f , ; > rk ; 5 ,V > ) ;-j .. r 3 

>. \r> c .,&■*!: 

)^tj] = d* ; s » . i .« v • ii; 0-. ».-;•* j I ;L. r-;k J 


The <d e t ermina t ion ^df A(t 0 ) J p , and t£ to satisfy the boundary ' condi- 
tions of" condition (d) subject ‘to "the initial Whditfohs (c6ndlti6h : (a)) 
and the differential equation constraints of cond’itiorf (c) is 1 tfie ISddttdiry 
value problem to be solved* Methods for its solution are discussed below* 

< > A Ky • ■» 1 i • . .... ; ... i i ,y ,-u > *< .. : » . i- *: i '• *r. - :■ ’ j ' •;.? 


-lJ •{■ 


6. Solution of Boundary Value Problem j 


■i ij l £i»'iO ,.'_i ; j 


As can be seen in equation., (107;),. the values of Mt 0 ) are the only 
variables that can be used to change the values of X(t) and A(t), because, 
by condition (a), X(t Q ) and t 0 are fixed, and, by condition (.c),, the 
differential equations and optimal u(t) relation must hold. Thus, the 
boundary conditions given in (d) of equation. (10j7), cap- be . written. as 

F[A(t 0 ),p,t f ] = 0, ' ' ' ' ' (108) 

, ^ ; . ; I, \ ' i ’>} . ( _ t i Z ^ 

where the dimension of F is equal to 'the dimens ion of G, plus the dimen- 
sion of A, plus one. Actually F has the form 


F[X(t f ),A(t f ),p,t f ] = 


G [X (t f ) , t f ] 


A (t f )+[0,0,0,0,0,0,0,l]+p 
p T G[X(t f ),t f ] : «•>:, 


Ti SG[X(t f ),t f ] 
dX(t f ) 


: : IV. i 


(109) 


56 



but X(t£) arid Mtjf) appearing in the right-hand side of the above equa- 1 
tion depend implicitly on A(t Q ) which leads to the functional form given > 
in equation (108). Now, for a guessed set of A(t Q ), p, and tf (denoted 
by A£,p*,t£), a Taylor series expansion can be written: 


+ . . . . 


* 

(HO) 


F[A(t o ),p,t f ] 


F[A*,p*,t*] + 


rdF[Mt o ),p,t f ]n 

BA(t o ),p,t f 


o’ p jt f 


AA(t o )j 

Ap 

At 


The higher order terms in the above expression can be neglected, a 
fudge factor K (where 0 < R i 1) can be chosen, and the above equation 
solved to yield the corrections (AA(t Q ), p, and Atf) which, when added 
to Ap, p* and tf and used to generate a new trajectory, will produce a 
decrease in the error term F* That: is, 


AA(t 0 ) 

Z\p 

Atf 


-K 


-BF[A(t o ),p,t f ] 
c3A (t Q ) , p , t^ 


-"1 


A jP ,t f 


F[A o ,p ,t f ] 


( 111 ) 


It can be seen that, as long as the matrix which must be inverted in 
the above equation is nonsingular, a value for K can be chosen as close 
to zero as necessary tp give a set of corrections (Z^A(t Q .) , Ap, and At^) 
which, when used to gene rate a new trajectory,, will make the square 
root of the error term (n/W F) less than or equal to 1 - K times the square 
root of the original error (/f t (Aq,p*, tg) F(A*,p*,t|) # Thus, if equa- 
tion (111) is used iteratively and the matrix invefse ; continues to 
exist, a set of values for A(t Q ) , p, and t^ can be found which will 
make F[A(t Q ) ,p, tf ] = 0. * , ... 


To obtain numerical values for the matrices 


-dF[A(t o ),p,t f ] i 


aA( to ) 


A O > P 9 ^ JT 


and 


r SF[A(t o ),p,t f ]-] 


-y* * . /V 3 

\’P .tf 
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the chain rule is used on the form of F given in equation (109) and 
results in the following equation: 


-dF[A(t Q ),p,t f ]- 


-dF[X(t f ),A(t f ),p,t f ]- 

/^ X (t f )\ 

» ( t o ) j 

X.* 

A„>P » t f 

dX(tf) J 

* W(t 0 >/* 


-dF[X(t f ),A(t f ),p,t f ]- 

dX(t f ) 


'dA(t f )^ 


( 112 ) 


~ciF[A(t o ),p,t f ] 


|-\ * Jc . * 

J A 0 »P 


F[A(t ) ,p, t f k* * 

u /x o » F ? 


(113) 


The matrix 


rSF[A(t Q ) ,p,t f ]-i 





f 


is obtained by explicit partial differentiation of the expression for 
F[A(t Q ) , p, tf] given by equation (109). Thus, the expression for 


-dF[A(t o ) ,p,t f ]- 
_ dA(t o ),p,t f 


A*, 


P*t 


* 

f 


can be written by combining the expressions in equations (112) and (113) 
with the expressions for 


rdF[A(t o ),p,t f ]- 

dp 


k * ■ k 

J A 0 ,P 
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In equation (112) the matrices 


■dF[X(t f ),A(t ),p,t ]-. rSF[X(t f ),A(t ),p,t ] 

. ^0“) and L SMV 


can be obtained by explicit partial differentiation, but the matrices 

^ x(t fb /dMO\ 

V* an W( C 0 V* 

must be obtained from the numerical integration of the following set 
of matrix differential equations: 



(114) 
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Note that condition (c) above allows 


3 “( t ) 

, , 3 A <t 0 ), . , n . ; SS . ? ,.K 

> .! : X * 

to be determined at each instant of time because the strengthened 
Legendre Clebsch condition requires the matrix 


d g X T f (X(t).u(t)) ! .1 5 ' ’■ >' r 
c)u(t) Su(t) : — 

to be negative definite. Thus, the above system can be integrated , 
numerically along with the system given in equation (107) , using the 
guessed values (A*, p*,tf) and the initial conditions 


Sx<t ) r ■ ; v ) dX(t ) ' ' _ ; ' ' ~ b 

ax(t o ) ~ t°i and sx(t o > " r hJ* ' i; >■ 

This completes the explanation of Newton’s method for the solution of 
the boundary value problem. 

When a solution of the boundary value problem by Newton’s method 
is obtained as discussed previous ly, 5 the sufficiency conditions given 
by conditions (b) jand (c) of r equation (95) can. be checked '.very easily. 
The normality condition will be satisfied if the matrix ^ /• 


SF[X(t o ),p,t f ]-|_ , 

_ :SX(t o ) ,P,t f j "j 

is positive definite 1 when evaluated for values of A(t Q ) ,p, and t^ which 
satisfy equation (107) . Reference 17 gives a very good explanation and 
a computer program listing for determining the definiteness of a matrix. 
The proof of this normality condition’s contribution to the sufficiency 
conditions requires a detailed and difficult examination of the second 
variitiotr as s-hoWh; in reference 16,; qnd wiif rip t; ;be' attempted (here., ) 


The easiest way to check a trajectory satisfying equation (107) for 
conjugate points is to generate the feedback guidance matrix. To do 

.J ' } ”• 3 ; ihZiHiXVXjl-': 1' ! + 
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this the differential equation for the transition matrices 


'dX(t) \ 


and 



are needed. These are obtained as equations (114) were obtained. 
That is. 


<a) 5T (Ho^) " ScfJ) ' sa " e as (lu) ,” ith 

replaced by X(t Q ) , 

(b) J£ (l^y) = liogr = same as (114) with 

replaced by X(t Q ), 

(c) [0] = same as (114) with A(t Q ) replaced b; 


The initial conditions for these equations are 

/dA(t )\ 

" 1 an,J Ur tt>) * t01 - 

Thus, for values of X(t Q ), p, and t£ which satisfy equations (107), the 
transition matrices 


3x(t o ) 

Sx(tJ 


Mt 0 ) 


A(t 0 ) 


(115) 


7 X(t ). 


^X(t) \ fbxj t) N 

^X(t o )/ W(t 0 )J 


bMt) \ 

^<t 0 V 


and 


a yJ 


where (t Q ^ t £ tf) can be obtained by numerically integrating equa- 
tions (114) and (115). Also needed are the vectors 


Ml l 


and 


M 1 

at o 
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which are defined by differential equations similar to equations (114) 
and (115). That is 


dX(t) 
' St o 

- same as (114) with 


replaced by t Q . 

SA(t) 
" at 

= same as (114) with 

o 

replaced by t Q . 


(116) 


(c) [0] = same as (114) with A(t Q ) replaced by t . 

Note that dX(t 0 )/dt 0 = X 0 and c)A(t 0 )/dt 0 = A Q . 

With all of these transition matrices, a set of trajectories about 
the trajectory-satisfying equation (107) can be represented by 


««=> - (S(?y) + M<t o> + At c 


(117) 


*<*> ' (fe) at(t o ) + (te) M(t o> + (^-) A V 


(118) 


Equations (117) and (118) state that small changes in the initial condi- 
tions which satisfy equations (107) are mapped by the transition 
matrices into small changes AX and aA at any time t where t Q S t § tf. 

In the set of trajectories defined by equations (117) and (118), a sub- 
set can be defined (all of which satisfy equation (107)) by requiring 
the AX(tf) and AA(t f ) produced by AX(t Q ), AA(t Q ), and At Q to satisfy 
condition (d) of equation (107). Thus, AF produced by substituting 
AX(tf) and AA(t^) into equation (109) must be zero. That is. 


AF = (dF/dX f ) AX f + (dF/dA f ) AA f + (dF/dp) Ap + F At f = 0. (119) 
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Substituting equations (117) and (118) (evaluated at tf) into (119) 
gives 



Equation (120) can be written as 


AA„ 

l—U 


r- -j 




Ap 

+ B 

At 

At 


o 

r 




where 


( 121 ) 


A = 








B 





Note that A is the same matrix which appears in equations (110) and 
(111) and must be positive definite from the normality condition. 
Thus, A is nonsingular so that 


AA 

o 


AX 

Ap 

P3 

H 

< 

If 

0 

At^ 

At f 


LT O J 


( 122 ) 
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Equation (122) can be partitioned to give 


where 


(a) 

AA 

o 

- Dl 

AX 

o 

+ 

d 2 

At: 

(b) 


= d 3 

AX 

o 

+ 

d 4 

At 

i 

(c) 

At ^ 

= D 5 

.AX 

o 

+ 

Ds 

At 


A" 1 B = 


Di 

d 5 


d 2 

d 4 

d 6 


(123) 


Now either Atf or At Q can be chosen to be zero. Usually, for purposes 
of comparison, it is better to choose Atf to be zero. Then, from equa- 
tion (123) part (c), it can be seen that 


At Q = -(1/D 6 ) D 5 AX o . (124) 


Equation (124) can be substituted into parts (a) and (b) of equation 
(123) to give 


(a) AA q = [Dj. - (1/D 6 )D 2 D 5 ]AX o = D T AX q 

(b) A*> = [D 3 - (1/D 6 )D 4 D 5 ]AK o = D s AX q . 


(125) 


Now part (a) of equation (125) and equation (124) can be substituted 
into equations (117). That is, 


AX(t) = 



AX(t o ) 



[d t ax q ] + 



[-(1/D 6 ) d 5 ax o ] 


(126) 


AA(t) 



«(t o ) 


•<$?) 


[D 7 AX ] + 

{ O 



[- (1/D 6 ) d 5 ax q ], 

(127) 
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Equations (126) and (127) can be rewritten to give 


AX(t) = P(t) AX(t Q ) 


(128) 

AA(t) = Q(t) AX(t Q ), 


(129) 

where 



'(*> - dSF) + 

x O y N o / 

” - 

(l/Ds) Ds 



(1/D 6 ) d 5 . 


The matrix P(t) which appears in equation (128) is the matrix which 
must be nonsingular for all t (where t Q ^ t < tf) in order to state that 
no conjugate points exist on the trajectory satisfying equation (107). 
Thus, this matrix can be computed at each time t and inverted if it is 
nonsingular (i.e., its determinant is nonzero). If there are no conjugate 
points on an optimal trajectory, then linear feedback guidance can be 
accomplished by solving equations (128) for AK Q and substituting the 
result into equation (129). That is, 


AA(t) = [Q (t) P _1 (t)] AX(t). 


(130) 


Also, the solution of equation (128) for AX(t Q ) must be substituted into 
equations (124) and part (b) of (125) to give 

Ap = [D s P _1 (t)] AK( t) (131) 

At Q = [- (1/D 6 ) D 5 P _1 (t)] AX(t) » (132) 


Then equations (130), (131), and (132) are the linear feedback guidance 
equations. That is, for a measured or given value of AX(t) (where 
t ^ t < t^), AA(t), Ap , and At 0 can be computed. Note that, with 
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values for AX(t) and AA(t) , equation (105) can be used to compute u(t) 
directly, or if equation (105) is linearized , then Au(t) can be computed, 
which is the usual form of feedback guidance. In the present MASCOT 
implementations, equations (130), (131), and (132) are not used, but 
they are now being examined for use in improving the computational 
speed of the atmospheric portion of MASCOT. 

The preceding discussion completes the solution of the boundary 
value problem by Newton's method and shows how the complete set of 
sufficient conditions can also be checked with little additional 
effort. As can be seen, this approach differs from the "sweep" method 
explained in chapter 6 of reference 16. After examining the computa- 
tional aspects and flexibility of both procedures, Newton's method 
(as explained here) has been selected as the most favorable approach 
for solving the boundary value problem and also for performing the 
sufficient condition Checks. More support for this viewpoint can be 
found in some of the statements and numerical results of reference 15. 

A considerable amount of effort is involved in establishing a com- 
puter program to perform a solution of the boundary value problem and 
make the sufficiency condition checks by Newton's method. Thus, 
simpler methods are usually used for problems which are in an early 
and formulative stage of development. One of the most effective 
examples of these simpler approaches for solving the boundary value 
problem is a modification of the secant method, which will be explained 
next. 


As before, the boundary values as given by equation (108) must be 
satisfied. To simplify the notation, a new vector will be defined as 


Mt 0 ) 


x = 


p 

fc f 


(133) 


Note that this X as defined above has no relation to any previous defini- 
tion of X. Then the solution of the boundary value problem is a vector 
X = (x!,x 2 , ..., x n ) that will satisfy equation (108) written in the 
following form: 


fl(Xl,X 2 , 

F(X) = f2 < x l’ X2 > 


x ) = 0 
x n ) - 0 


• • • , 


( 134 ) 



In the explanation to follow, for a given set of values of the X' s , 
the f f s can be evaluated only by numerically integrating the differential 
equations with their given initial conditions and the given set of values 
of the X f s completing the initial conditions. If a coasting phase occurs 
at any point during this numerical integration procedure, very rapid 
closed-form solutions are available {8, 9, 10, 11] to obtain the trajec- 
tory across this coasting arc without the need for time-consuming 
numerical integration. Then after the integration is complete, the f*s 
can be evaluated. There are many techniques available for numerically 
integrating ordinary nonlinear differential equations (such as Runge- 
Kutta and Shanks formulas), but the most effective yet developed are 
given in references 18 and 19. Thus, with a given numerical integra- 
tion technique, the functions f can be evaluated for each set of X* s . 

As can be seen in the explanation to follow, the transition matrix dif- 
ferential equations (113) are not needed for this approach, which accounts 
for most of its simplicity. 

This method begins with a reasonable but nonetheless arbitrary first 
choice of x l9 x 2 , ..., x n . A variation in each of the x T s must be made. 

One possibility is to define 


y. =k.x., i = 1, ..., n and 0 ^ k. < 1 

h i x 9 ’ i 


and then set up the matrix [X] : 


(135) 


” — 


— — 

*1 • • • x n 


x n x i2 ... x in 

Y 1 x 2 ... x n 


X 21 x 22 • • • x 2n 

Xi y 2 ... x n 

• 

• 

= 

• 

9 

9 

• 

• 

x i x 2 ... y n 


* 

X jl X j2 * * ' x j n 


(136) 


Note that j =* n + 1, and hence [X] is an n 4- 1 X n matrix. 

Now, for each of the j rows of [X], X^, evaluate the n f 1 s . This 
results in a j X n matrix of f values. Set 


fl(Xi) 

f 2 (Xi) 

... f n (Xx) 


f ll 

f 12 ... 

fin 

fl(Xg) 

• 

f 2 (X 2 ) 

... f n (X 2 ) 


f 21 

f 22 

f 2 n 

• 

ji(Xj) 

f 2 (X j ) 

... f n (xp 


f ji 

fj 2 ... 

f jn 


( 137 ) 
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In order to make [X] square (for reasons to be seen), make U the 


[X] = 


[X] 


u 


Ls a 

column 

vector 

of 

j I's, 

and define 

*ll 

*12 - 

. X ln 

1 



*21 


• x 211 

1 



• 

• 



• 

(138) 

X jl 

X j2 •• 

* X jn 

1 




The next step in the procedure is to try to replace the row of [X] 
which resulted in the "worst" row of [F] (for example, the one having 
the maximum sum of the squares of the fj-, assuming proper normalization) 
by an X new which results in an improved F. To do this, first note 
that a matrix [A] exists which has the property 


[X] [A] = [F], 


(139) 


and set 


X 


new i 


1 


[A] 


= F 


new 5 


(140) 


where X new is a vector of the unknowns that should produce F new and 


F * C 
new 


old : 


0 ^ C < 1 


(141) 


and F old is a row of [F] selected by virtue of its having the minimum 
sum of squares of the f^. Note here that an X new is desired which will 
lead to the "worst" row of F being replaced by one better than the "best" 
row of [F] „ 

From (139), 


[A] = [X ]- 1 [F], 


(142) 
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and therefore (140) becomes 



new i 


i 


[X]- 1 [F] = F . 

new 


(143) 


The solution for X„^ TT proceeds as follows: Let 

new L 


P 


T 


X 


i 

i 

i 


new i 


i 


1 


[X]- 1 , 


(144) 


where P is a j = n + 1 column vector. Then, using (144) in (143), 


T mm 

P [F] = F or [F] P = F 
L new new 


(145) 


Equations (145) are a set of n equations in j = n + 1 unknowns. However, 
from (144) 


T 

P [X] 



1 "1 


rX -| 

= 

x ! l 

or [X] T P = 

_new_ 


new i 

L l J 

1 


(146) 


T 

Looking back at (138), note that [X] has a bottom row of l T s, so 
that from (146), the following equation is obtained: 


j 

^ - 1. (147) 

i=l 
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This provides the jth equation needed to solve for the j unknowns 
above, as a system of simultaneous equations; i.e.. 


[F] T P = F T and > P. = 1 (148) 

new i 

i=l 


are sufficient to obtain the j components of P. 
This may also be written 



Pi 



r ^T i 
F 


p 2 


T 

r F i 
new 

“"t" 

L u J 


• 

• 


_ 1 - 


p. 

j 




rp 

where U is a row of vector of j l r s. 


(149) 


The nonexistence of the inverse in (149) does not necessarily imply 
that the solution as a set of simultaneous equations does not exist, 
since the rank of the augmented matrix may be the same as the original 
one # 

Once P is known, (147) gives X new , since from (146) 


T 

[X] P = 


-x 

_new 

L 1 J 


(150) 


This completes the discussion of the methods for solving the 
boundary value problem as given by equation (107). Several computer 
programs are available from the authors for solving the boundary value 
problem as given by equation (107). The listing of these computer 
programs are not included here because they are almost constantly 
being modified. Some preliminary results on space shuttle type 
vehicles have been obtained with the secant method (as explained here) 
and they are discussed in the Results section of this report. In refer- 
ence 9, a computer program listing for a Newton 1 s method solution of 
the Vacuum Trajectory Optimization problem is given. This program, 
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called SWITCH, is available for use. In reference 8 some examples of 
the application of feedback guidance associated with Newton 1 s method 
are given to show that feedback guidance matrices can be used for 
guidance and to check for the satisfaction of sufficiency conditions. 


D. Conclusion of Mathematical Development 

The mathematical development presented here for the- MASCOT concept 
has of necessity been rather complicated. It is hoped that the reader is 
able to spend enough time studying the ideas presented so that he can 
understand them. This understanding of the general idea and the numeri- 
cal results presented in the next section should convince most readers 
that the MASCOT concept as outlined here can be developed into an effic- 
ient, effective, and flexible guidance scheme for vehicles of the space 
shuttle type. 

To reiterate and emphasize the key mathematical ideas discussed in 
this section which now make the onboard solution of the trajectory 
optimization problem (MASCOT) an attractive guidance scheme for space 
shuttle type vehicles, the following statements are listed: 

First, the unified set of guidance equations and performance 
criteria will reduce the computer storage requirements for onboard 
implementation. 

Second, the techniques for obtaining a solution of the boundary- 
value problem have been carefully studied and the so-called "shooting 
method" selected because of its speed and reliability. The shooting 
method algorithm is easy to program and is compact in size. The 
alternatives to the shooting method, such as steepest descent and quasi- 
linearization, must use stored time functions so that the size and com- 
plexity of the algorithms are increased. 

Third, the development of the Fehlberg-type, Runge-Kutta numerical 
integration routines has significantly speeded up the numerical integra- 
tion process. Further improvements to numerical integration techniques, 
such as the variation of parameters method of reference 20, can also be 
easily incorporated into the technique. 

Fourth, since the boundary value problem can be considered as the 
problem of finding. the solution to a system of simultaneous nonlinear 
equations, any improvements in these techniques which may occur in the 
future can again be incorporated easily into the technique. 

Finally, the many improvements in computational hardware which 
have resulted in flight-worthy high-speed digital computers also adds 
to the attractiveness of the concept. 
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VII. NUMERICAL RESULTS 


The resources for developing the OPGUID, SWITCH, and MASCOT algo- 
rithms have been limited* To obtain the maximum benefit from available 
funds and manpower, the effort has been concentrated in the area of 
problem formulation and improving the numerical techniques needed to 
solve the resulting boundary value problem rather than in performing 
detailed guidance simulations. Therefore, only limited results are 
available at this time. These will be discussed in the following 
paragraphs. The results of the detailed guidance simulation studies 
now being performed will be published at a later date. 


A* OPGUID 


The feasibility of solving the trajectory optimization problem in 
real time for onboard guidance was demonstrated in 1966 by Brown and 
Johnson [7] for an ascent to orbit mission (vacuum). Since these 
results can be found in t he literature, they will only be summarized 
here. 


The mission selected was ascent to a circular orbit with the S-IVB 
stage of the AS-204 vehicle. To evaluate the stability and optimality 
characteristics of the real-time optimal guidance scheme (OPGUID) under 
realistic conditions, a number of simulations were made and compared to 
results obtained, making the same simulations with the IGM equations. 
These included thrust and mass rate variations of ±5 percent and engine 
mixture ratio shift times of ±30 seconds from nominal. OPGUID used a 
variable cycle time ranging from 10 seconds initially to one second 
near the end. In all cases, the OPGUID scheme delivered more payload 
to orbit than the IGM scheme; however, these differences in all cases 
were small (2 to 70 kg), and the main point to be made from this 
study was that it demonstrated the feasibility of using OPGUID as an 
onboard guidance scheme. The computer speed and storage requirements 
for OPGUID can be met by the Saturn Launch Vehicle Digital Computer; 
thus, OPGUID could serve as an onboard guidance scheme for the Saturn V 
vehicle. 


B. SWITCH 


Recently, under NASA contract NAS8-21315, a sophisticated multi- 
burn optimization program, SWITCH [9], was developed which appears to 
have the speed and reliability needed for onboard guidance signal 
generations. The convergence properties of the SWITCH program were 
demonstrated by applying it to a multi-burn rendezvous mission. The 
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target vehicle was located in a synchronous orbit with an altitude of 
19,300 n.mi., and the pursuit vehicle was assumed to be in a 100 n.mi. 
circular orbit. The number of separate burn arcs was set at two 
(coast-burn-coast-burn) and the initial values of the co-state vector 
and the lengths of arcs were chosen to correspond to values obtained in 
a previous study. The initial values of independent variables were 
extremely close to the required values, and the SWITCH algorithm con- 
verged to the solution in one iteration. Since the initialization was 
essentially equal to the converged solutions, the convergence properties 
were evaluated by deforming the final orbit into an ellipse with the 
eccentricity magnitude varying from 0 to .5. The pursuit vehicle orbit, 
vehicle parameters, and initial estimates of co-state and arc times 
were fixed. The results are taken from reference 9 and are shown in 
Table 1. Notice that the worst case required only five iterations. 


Table 1. Low Altitude to Coplanar Synchronous Rendezvous Mission 


Target Orbit 
Eccentricity 

1st Coast 
Arc 
(sec) 

1st Burn 
Arc 
(sec) 

2nd Coast 
Arc 
(sec) 

2nd Burn 
Arc 
(sec) 

Number of 
Iterations to 
Converge 

0 

399.83 

255.82 

18,729.51 

129.11 

i 

.05 

401.16 

253.47 

17,572.70 

134.87 

3 

.1 

402.48 

251.14 

16,541.83 

140.57 

3 

.5 

412.62 

232.95 

11,155.52 

184.04 

5 

First itera- 
tion speci- 
fications 
for all 
cases 

400 

255.82 

18,727.46 

1 

129.11 ! 



Although the SWITCH algorithm has not been subjected to the stringent 
guidance simulations that are necessary to prove the flight worthiness of 
an onboard guidance scheme, it is felt that the excellent convergence 
properties of the SWITCH algorithm and the speed per Iteration (one quarter 
second per iteration on a CDC 6600 computer) of the algorithm make it 
applicable for onboard guidance. The computer requirements for SWITCH 
are greater in terms of speed and storage than for OPGUID; however, they 
are not considered to be excessive. It is estimated that the SWITCH 
algorithm would require 7,000 storage locations for onboard implementation. 
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Assuming a state-of-the-art type of flight computer, the guidance cycle 
time (time for complete converged computation of new optimal guidance 
command time history) would be approximately 2 seconds. 


C. MASCOT 


The MASCOT guidance algorithm is an extension of the SWITCH and 
OPGUID formulations. It can be reduced to the SWITCH algorithm by elim- 
inating the aerodynamic forces, and it can be reduced to the OPGUID 
algorithm by eliminating both the aerodynamic forces and the multi-burn 
logic. Therefore, the results presented for OPGUID for ascent to orbit 
(vacuum), and the orbit transfer and rendezvous results obtained with 
SWITCH are indicative of MASCOT for these types of missions. Results 
for atmospheric flight are preliminary and it remains to be shown that 
the MASCOT scheme possesses the speed and reliability for onboard guidance 
for atmospheric flight. As mentioned earlier, studies to determine the 
properties of MASCOT are underway and will be documented as soon as results 
are available. 

The need for a shuttle guidance scheme with the sophistication of 
MASCOT is demonstrated by the performance results presented in Table 2. 

TABLE 2 


Case 

Time of 
Max Q 
(sec) 

Max Q 
(kg/m 2 ) 

Angle of 
Attack 
at Max Q 
(deg) 

Cutoff 

Weight 

(kg) 

A Cutoff 
Weight 
(kg) 

A 

IGM from 
Staging* 

63 

2310 

0 

122,259 

- 

B 

Optimal from 
12 sec. 

46 

4087 

8.0 

123,380 

1121 

C 

Optimal from 
70 sec. 

63 

2310 

0 

122,519 

260 

D 

IGM from 
70 sec. 

63 

2310 

0 

119,706 

-2553 


All trajectories are zero angle of attack until initiation 
of IGM on optimal shaping. 
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The data shown in Table 2 indicate the effects on payload of optimally 
shaping the space shuttle trajectory during the atmospheric region of 
flight. Case A represents trajectory and vehicle data for a shuttle 
vehicle that was flown through the atmospheric portion of flight with 
a zero angle of attack, and the IGM steering logic was used to shape 
the remaining portion of the trajectory. Case B represents data 
obtained by shaping the trajectory with the MASCOT logic from 12 sec 
after lift-off through orbit insertion. This gain in payload of 1121 kg 
may not be fully obtainable because shaping the trajectory during the 
atmospheric portion of flight will result in increased structural loads. 
Case C represents data that were obtained by shaping the trajectory with 
the MASCOT logic from 70 seconds (7 sec after maximum dynamic pressure) 
after lift-off through orbit insertion. This increase in cutoff weight 
of 260 kg is potential payload since the maximum dynamic pressure is 
the same as in case A. More important than the 260 kg increase in pay- 
load is the fact that using MASCOT to guide the vehicle from 70 sec 
(or earlier) provides closed-loop guidance early in flight, and thus 
the effects on payload of perturbations (winds, thrust variations, etc.) 
during the first 70 seconds of flight can be minimized. The results of 
using IGM for closed-loop guidance at 70 sec is shown in Table 2 as 
case D. This loss in payload of 2553 kg is due to the following: (1) 
atmospheric forces are not modeled in IGM, and (2) the optimal steering 
angle during the atmospheric portion of flight is nonlinear (IGM is 
a linear steering law). Shown in figure 12 is the optimal steering 
angle versus time for case B. Notice the nonlinearity from 12 sec to 
approximately 130 sec. This indicates that linear steering laws are 
far from optimal and should not be used for guidance in this region of 
flight even if atmospheric forces are included in the guidance equations. 

There are no significant results to report for the reentry phase. 
Early computations have shown a very high sensitivity of the reentry 
trajectories to the choice of the initial values of the Lagrange multi- 
pliers and consequently a very poor convergence rate. However, recently 
the approach outlined on pages 25 and 26 has been implemented in order 
to eliminate this high sensitivity and poor convergence. This method 
assumes that u(t) could adequately be represented as a linear function 
of time during the atmospheric portions of flight. 

Studies are currently underway to determine whether it will be neces- 
sary to use the solutions so obtained as first approximations to the true 
optimum. Presently, no further corrections appear to be necessary. 
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Figure 12. 


VIII. CONCLUSIONS 


Fast, efficient, compact trajectory algorithms for a total trajec- 
tory from lift-off to landing have been developed using a numerical 
solution of the optimization boundary value problem. These algorithms 
use a rapid numerical integration routine and an efficient iterative 
solution of the simultaneous nonlinear equations involved. The imple- 
mentation of one of these algorithms (using the available vastly 
improved on-board flight computers) would accomplish a major advance 
in guidance techniques. That is, preflight analysis would be reduced 
to a minimum, thus making it possible to achieve Shuttle launches with 
airline type of turnaround times. Also, special tuning constants and 
functions (which have heretofore been necessary to account for approxi- 
mations introduced to simplify onboard computations ) would be eliminated. 
Thus, this proposed guidance system would produce (on board the vehicle) 
optimal trajectories for an extremely wide variety of missions and 
vehicle characteristics. 
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APPENDIX I 

COMPUTER PROGRAM LISTING 
FOR PARAMETER OPTIMIZATION DECK 
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262 100 RXVQB CI 3 *RXV 0 B CI 3 /RXVO 

263 RVR [ 13 «RXV 0 B [ 2 ] *R 09 [ 33 *RXV 0 B [ 3 ] *R 0 B [23 

26 ^ RyR C 23 =RXV 0 B [33 *R 0 B C 13 -RXVOB [13 *R 0 B C 33 

265 RVR [33 =RXV 0 B [ 13 *R 0 B [ 23 *RXV 0 B [ 23 *R 0 B [13 

266 C START OF TRIAL INTEGRATION 



267 T 8L a *05 

268 CALL DEQ [X,TI,DX] 

269 PRINT 260 

270 CALL PRINT CTIjX#DX3 

271 TsTI 
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